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ABSTRACT 

We investigate how star-forming galaxies typically assemble their masses at high redshift. Taking 
advantage of the deep multi-wavelength coverage of the GOODS data set, we select two of the largest 
samples of high-redshift star-forming galaxies based on their UV colors and measure stellar mass of 
individual galaxies. We use template-fitting photometry to obtain optimal estimates of the fluxes 
in lower-resolution ground-based and Spitzer images using prior information about galaxy positions, 
shapes, and orientations. By combining the data and realistic simulations to understand measurement 
errors and biases, we make a statistically robust determination of stellar mass function (SMF) of the 
UV-selected star- forming galaxies at z ~ 4 and 5. We report a broad correlation between stellar 
mass and UV luminosity, such that more UV-luminous galaxies are, on average, also more massive. 
However, we show that the correlation has a substantial intrinsic scatter, particularly for UV-faint 
galaxies, evidenced by the fact there is a non-negligible number of UV-faint but massive galaxies. 
Furthermore, we find that the low-mass end of the SMF does not rise as steeply as the UV luminosity 
function («uvlf ~ — (1-7 — 1.8) while «smf ~ —(1.3 — 1-4)) of the same galaxies. In a smooth and 
continuous formation scenario where star formation rates are sustained at the observed rates for a 
long time, these galaxies would have accumulated more stellar mass (by a factor of ~ 3) than observed 
and therefore the stellar mass function would mirror more closely that of the UV luminosity function. 
The relatively shallow slope of the SMF is due to the fact that many of the UV-selected galaxies are 
not massive enough, and therefore are too faint in their rest-frame optical bands, to be detected in 
the current observations. Our results favor a more episodic formation history in which star formation 
rates of low-mass galaxies vary significantly over cosmic time, a scenario currently favored by galaxy 
clustering. Our findings for the UV-faint galaxies at high redshift are in contrast with previous studies 
on more UV-luminous galaxies, which exhibit a tighter SFR-M sta r correlation. The discrepancy may 
suggest that galaxies at different luminosities may have different evolutionary paths. Such a scenario 
presents a nontrivial test to theoretical models of galaxy formation. 

Subject headings: galaxies: evolution — galaxies: formation — galaxies: high-redshift — galaxies: 
starburst - surveys 



1. INTRODUCTION 

Our knowledge of the distant universe has expanded 
substantially over the last decade thanks to deep 
multi-wavelength surveys. Currently, galaxies are 
routinely identified out to z ~ 7 and beyond, and 
their key quantities such as star formation rates (SFR), 
stellar masses (Mst ar ), and extinction are measured via 
spectroscopy and multi-wavelength phot o metry (e.g, 
Papovich etail I200H IShaplev et all 120011: iStark et al.l 
2009t iMannucci et al.1 12009() . Large- area surveys have 



also enabled us to characterize the statistical prop- 
erties of the high-redshift galaxy population as a 
whole. The UV luminosity function primarily reflects 
the distribution of star formation rates within the 
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Gonzalez et al.ll2011[) . Morphological studies have quan- 



tified the rel ative spatial distribut i on and size of stellar 
components dFerguson et all 120041: iBouwens et al.l [20041 : 
iRavindranath et all 120061: lLotz et al.l I2006D . Finally, 
the galaxy clustering measures the large-scale spatial 
correlation between galaxies and how it compares with 
that of underlying dark matter and dark matter halos 
(iGiavalisco fc Dickinson! 120011: lOuchi et al.ll2004bl 120051: 
Adelberger et alj 120051 iQuadri et al.l 120071: iLee et atl 
2006L I2009af l7~ Careful comparison of these statistical 



quantities at different cosmic epochs can provide power- 
ful tools to gain insight into the evolutionary sequence 
of galaxy assembly in the high-redshift universe. 

This tremendous observational progress offers the op- 
portunity to establish whether galaxy growth is dom- 
inated by smooth accretion or more episodic events. 
While both modes certainly exist and are observed 
at high redshift, the relative importance of the two 
modes is currently unconstrained at z > 3 (but see 
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Rodighicroiril] [20111: lElbaz et al.l [201ll for measure- 
ments out to z ~ 2). Furthermore, there are observa- 
tional evidence that the dominant mode may va ry with 
galaxy or halo mass (e.g. , Ivan der Wei et aLll2011f ) . A key 
measurement that directly constrains the main mode of 
star formation, or average star formation history, is the 
location of galaxies on the SFR-M sta r plane. If galaxies 
assemble smoothly over a time scale comparable to the 
Hubble time, the majority of galaxies woul d form a tight 
seque nce with a slope of close to unity (Noes ke et al.l 
2007). On the other hand, in the case of bursty/episodic 
star formation, one should observe significant scatter 
about the mean relation. 

Unfortunately, at high redshift, different studies report 
conflicting resu lts. Some find a strong correlation with 
a tight scatter (IDaddi et al.ll2007U IPannella et al.ll200!i 
IMairdis et all I2010t ILee et al.l 1201 11 ISawickil 1201 lh while 
others find no or weak correlation with larger scatter 
(Sha plev et~aTll2005t Heddv et al.|[200l iMannucci et al.l 
2009). Some of this discrepancy may be real, if the main 
mode of ga laxy growth changes with stellar mass or host 
halo mass (jRenzinil 12001 ILee et al.|[2009al ). or with dif- 
ferent "types" of galaxies. However, it is likely that the 
discord between different measurements may be at least 
in part due to different and poorly understood selection 
effects suffered by different samples and/or different cal- 
ibration methods employed by these studies to derive 
various star formation rates and stellar mass. 

Some of the existing samples of galaxies are based on 
a UV color-selection while some others use a combina- 
tion of UV and optical colors. Some even require an ad- 
ditional detection in the mid-infrared (such as Spitzer) 
complicating the interpretations of their results in the 
global picture. Calibration is another issue at hand. For 
example, SFRs derived from polycyclic aromatic hydro- 
carbon (PAH) emission are known to depend strongly 
on galaxy metallici ty and the interna l distribution of HII 
regions and dust (|Calzetti et al.ll2007ft even when the pre- 
cise redshift is known. In addition, the presence of active 
galactic nuclei (AGN) in some of high-redshift galaxies 
may contribute to the ob served flux in the mid- infrared 
(e.g., IDaddi et al.ll2007al) . and subsequently add to the 
observed scatter in the SFR-M sta r relation unless such 
galaxies are identified or the AGN contribution is ac- 
counted for. These complications make it challenging 
to constrain the intrinsic SFR-Af star scaling law even for 
galaxies with the most comprehensive multiwavelength 
and spectroscopic coverage. Furthermore, at z > 3, even 
these diagnostics become out of reach as the PAH emis- 
sion redshifts to the longward of 24 um (and Ha redshifts 
to > 3/im; but see Shi m et al.ll201lD . and only galaxies 
with the most extreme starbursts such as sub-mm galax- 
ies can be detected via other Spitzer MIPS bands or Her- 
schel. 

Here we take an alternative approach to constrain the 
main mode of star formation by using two statistical 
measurements of galaxies, namely, the UV luminosity 
function (UVLF) and stellar mass function (SMF) of the 
same galaxies. The former constrains the distribution of 
recent star formation in galaxies, while the latter mea- 
sures the distribution of existing stellar mass in the same 
galaxies directly tracing the integral of all past SF activi- 
ties. Hence, these two quantities represent the statistical 
distribution of the SFR (modulo dust) and stellar mass 



rather than those in individual galaxies. In essence, com- 
parison of the two distributions will shed light on how one 
traces the other (when the effect of dust extinction is ac- 
counted for). If they trace each other closely, the overall 
shape of the two distributions will be similar. Alterna- 
tively, any difference between the two distributions can 
be attributed to the statistical mapping between the two 
as the two quantities are measured from the same galaxy 
population. 

There are several advantages of this approach over 
the direct method discussed previously. First, as 
we will show later, the two quantities, UVLF and 
SMF, can be determined very robustly by combining 
the observations and extensive simulations via rigor- 
ous statistical analyses. Second, our approach bypasses 
the need to measure star formation rates very accu- 
rately for individual galaxies. As mentioned previ- 
ously, the estimation of star formation rates is hard 
even when precise photometry and redshift informa- 
tion are available, and has been the primary challenge 
to make progress. Instead, we take advantage of the 
deepest surve ys to statistically map the UV luminos- 
ity to SFRs (IBouwens et~ail 120091: IReddv et al.l IliiH: 
iCastellano et al.l 12011b iFinkelstein et al.l 120111 ). In par- 
ticular, at z > 3, this approach offers new possibilities to 
obtain useful insights into the typical mode of galaxy 
growth as deep imaging surveys are being conducted 
in the near- and mid-infrared from space (e.g., Spitzer 
Extended Deep Survey, Cosmic Asse mbly Near-infrared 
Deep Extragalactic Le gacy Survey; iGrogin et al.l 1201 lb 
iKoekemoer et al.ll2011h . 

In this paper, we investigate how star-forming galax- 
ies assemble their mass in the first two billion years. 
We take advantage of two of the deepest survey fields 
with the most comprehensive multi- wavelength coverage, 
namely, GOODS-N and GOODS-S field. The combi- 
nation of the large areal coverage (~0.1 degree 2 total) 
and depth achieved in these fields allows us to identify a 
large number of high-redshift star-forming galaxies and 
robustly measure their population properties from inte- 
grated galaxy light. We use the TFIT template-fitting 
photometry package to carry out the photometry and 
present the most extensive set of simulations to date 
to validate the techinque and quantify the uncertain- 
ties. As we will demonstrate later, TFIT greatly im- 
proves the photometric accuracy for individual measure- 
ments over the conventional methods as well as dramati- 
cally increases the number of sources with reliable multi- 
wavelength photometrjQ. 

Throughout this work, we use (f2 m , Oa, cs; ^ioo) — 
(0.28 , 0.72, 0.9, 0.72). M agnitudes are given in AB sys- 
tem (Okc fc Gunn|[l983T ) unless noted otherwise. 

2. DATA, SAMPLE SELECTION, AND 
MULTI- WAVELENGTH PHOTOMETRY 

2.1. Data 

The data set we use for our analyses consists of those 
obtained as part of t he Great Observatories O rigins Deep 
Survey (GOODS; IGiavalisco et~aLl l2004bD . The two 

9 For example, IStark et all 120091) noted that only 35% of the 
sources at z ~ 4 are sufficiently isolated to perform robust aperture 
photometry. On the other hand, this work using TFIT rejects less 
than 5% due to blending. 
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fields, GOODS-N and GOODS-S, each covers roughly 
160 arcmin 2 , include deep HST/ACS S435V606^775^850 
(F435W, F606W, F775W, and F850LP) and deep Spitzer 
IRAC [3.6/H, [4.5/im], [5.8/zmJ3 and MIPS [24/xm] 
data. The GOODS-S field is covered by deep J, H, 
and -ftTg-band data taken with the ISAAC camera, and 
the L/-band taken with the Visible Multiobject Spec- 
trograph (VIMOS) o n the Very Large T elescope (VLT: 
iRetzlaff et all 120101 iNonino et all 120091 respectively). 
The GOODS-N field is covered by the J and K s - 
band data taken with the Wide-Field Infrared Cam- 
era (WI RCAM) on the Can ada-France-Hawaii Telescope 
(CFHT; iWang et al.l |2010L L. Lin et al., in prep), and 
the t/-b and data taken wit h the MOSAIC prime focus 
camera (|Capak et al.l [20041 on the Kitt Peak National 
Observatory 4-m telescope. 

The 5cr limiting magnitudes of the data are tabulated 
in Table [TJ Therein, the limiting magnitudes for the ACS 
bands are computed within an 0'.'2 diameter aperture. 
For the rest of the photometric bands, we compute the 
photometric uncertainties (la errors) within the isopho- 
tal apertures defined in the detection band (z 8 5o-band; 
see £12.41 for details). 

2.2. "Dropout" Samples at z ~ 4 and 5 

We adopt the Lyman-break technique (iSteidel et al.l 
1996) to identify high-redshift star-forming galaxies 
("Lyman Break Galaxies" or LBGs, hereafter) at z ~ 4 
and ~ 5 with the colo r selec tion criter i a discussed i n 
IGiavalisco et~aT1 (|2004al) and iLee et al.l (|200l l2009al ). 
Extensive spectroscopy campaigns have shown that these 
selection methods are very robust with minimal c ontami- 
nation by interlopers (|Vanzella et al.ll2006l l2009f ) . While 
the adopted selection criteria differ slightly from o thers 
found in the literature fe.g.. iBouwens et al1l2007t ). our 
main conclusions should remain identical to those de- 
rived using different selection as the proper account of 
the selection efficiency will make appropriate corrections 
to compensate for the difference . We us ed the UV color 
criteria used by IGiavalisco et al.1 (|2004a| ) to select galax- 
ies at z ~ 4 and z ~ 5 as 

(B435 - Veoe) > 1.2 + 1.4 x (V 6Q6 - z 850 ) A 
(S435 - Veoe) > 1-2 A V 6 Q 6 - z S50 < 1.2, 

and 

(^606 ~ «775) > 1-5 + 0.9 X (i 775 - Z 850 ) V 
(VaoG - i 775 ) > 2.0 A (V &m ~ i 77b ) > 1.2 A 

(«775 - ^850) < 1-3 A 

S/N(B 435 )<2, 

where the symbols V and A are the logical "OR" and 
"AND" operators, respectively. In both samples, only 
sources with S/N(z 850 ) > 6 are considered. In addi- 
tion to the color criteria, we also removed sources that 
have stellarity inde x (using the CLASS_ST AR parameter 
in the SExtractor; iBertin fc Arnoutslll996l ) greater than 
0.8 from the sample when the source is brighter than 
26.2 mag in the detection band (z 850 -band). For sources 
fainter than this limit, our simulation suggests that the 
stellarity measurement is not as reliable. 

10 We note that the majority of our samples are not detected in 
the [5.8^tm] and [8.0^tm] bands 



The number of photometric candidates satisfying these 
selection criteria is 3088 and 987 at z ~ 4 and 5, re- 
spectively. We refer to these samples as £?435-band 
and V606-band dropouts hereafter as their targeted red- 
shift range require them to "drop out" in those pass- 
bands. Finally, we removed the most likely low- 
redshift inter l opers based on photometric redshifts from 
iDahlen et all (|20Toh . which are calibrated against the 
largest compilation of spectroscopic redshifts available in 
the GOODS fields (IWirth et al.l l200l ICowie et al.ll200l 
iVanzella et al.ll2006l I2009L Stern et al., in prep). Using 
their estimates, we remove sources whose total integrated 
probability P(z < 3) is greater than 70%. The number 
of sources removed using this criterion is 136 and 131 
ga laxies at z ~ 4 and 5, respectively. Figure 15 and 18 
of IDahlen et al.1 (|2010l ) show the distribution of photo- 
metric redshifts for the -B435— and Vgo6 _ band dropout 
samples. 

2.3. Spitzer MIPS -detected Sources 

We identify sources with the Spitzer MIPS 24 /im de- 
tection (with the formal signal-to- noise ratio S/N> 3). A 
bona-fide high-redshift galaxy at z > 3 with 24 /im de- 
tection (sampling the rest-frame A w 4 — 6 /im) implies 
the presence of unusually strong hot dust components 
produced by either an AGN or an extreme starburst 
(SFR> several hundreds M Q yr _1 ). Indeed, of the -B435- 
and V(506-band dropouts that are genuinely detected at 
24 /im, several are detected in the deep X-ray data 
([Alexander et a l. 2003; Luo et all2008|) , further support- 
ing the possibility that they harbor an AGN. We cross- 
matched the GOODS MIPS catalogs with the dropout 
lists using a tolerance of 0'.'5 to minimize blending prob- 
lems. For the £?435-band dropouts in the GOODS-N, we 
find 24 sources with MIPS detections (four of which are 
detected in X-ray). Of those, 11 objects have spectro- 
scopic redshifts (3 are foreground interlopers at z < 3 
and the other 8 at 3.403 < z < 4.604). In the GOODS- 
S, we find 24 MIPS detection (five of which are X-ray 
detected). Of the total of 7 objects with spectroscopic 
redshifts, one lies at z = 2.797 while the remaining 6 
lie at 3.055 < z < 3.891 on the low side of the redshift 
distribution for the f?435-band dropouts. Four of the six 
sources at 3.055 < z < 3.891 have X-ray detections. In 
summary, 44 out of 48 MIPS detections in the combined 
GOODS-North and South field have f(24 /im) < 82 /xJy. 
Three spectroscopically confirmed i?435-band dropouts 
have 100 /xJy< /(24 /im) <160 /iJy and are also X-ray 
sources, thus likely AGN. 

As for the V606 _ band dropouts in the GOODS-N, we 
find 4 MIPS detections none of which has spectroscopic 
redshifts available. However, two of them are very bright 
in 24 /im (113 /xJy and 347 /iJy), and thus are very 
likely to be foreground interlopers. In the GOODS- 
S, we find 11 MIPS detections three of which are con- 
firmed to be foreground interlopers (2=1.324, 1.981, and 
3.513) while another lies at z = 4.762 on the low-side 
of the redshift distribution for the V"606 _ band dropouts. 
The gal axy at z = 4.762 is also a known submillimeter 
galaxy (Cop pin et al.l 12009). Of those without spectro- 
scopic redshifts, one is unusually bright (109 fxJy) to be 
at z > 4.5. In summary, for the combined GOODS-N 
and GOODS-S y 606 -band dropouts, 12 out of 15 MIPS 
detections have /(24 /im) < 60 /iJy. Two with no spec- 
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TABLE 1 

The 5cr Limiting Magnitudes of the Data Set 
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Note. — For the ACS bands, the limiti ng m agnitudes are computed within an 0'.'2 diameter aperture. For the rest, isophotal aperture 
defined in the detection band is used (see i|2,4l l 



troscopic confirmation have /(24 /im) > 110 /zJy (no 
X-ray detection) and another has /(24 fj,m) = 347 /xJy 
(also without spec-z or X-ray). 

Based on these statistics, we conclude that rather few 
bona fide 5435-band dropouts, and almost no V(306-band 
dropouts are detected at 24 /zm. Thus, we exclude all the 
sources with MIPS detection from our sample. However, 
we note that inclusion or exclusion of these sources makes 
little difference in our main conclusions as their number 
is few. 

2.4. Template Fitting Photometry (TFIT) 

Reliable estimates of physical parameters begin with 
reliable photometry. This is a challenge for mixed- 
resolution datasets like GOODS, where the PSF width 
varies by an order of magnitude across the full wave- 
length range. Galaxies that are cleanly identified in the 
HST images may be horribly blended with their neigh- 
bors in the Spitzer images. 

Template Fitting (TFIT) photometry is designed to 
overcome these challenges and provide accurate flux 
estimates and colors for galaxies in multi-wavelength 
mixed-resolution data set such as the data used 
in this work. Whil e w e refer interested r eaders to 
iPapovich etaD (|20Qll ) and lLaidler et all (|2007l ) for more 
details of the TFIT algorithm, we also note that similar 
algorithms have been impl emented by others to obtain 
optimal photomet r y Ce.g.. iFernandez-Soto et al.l IT999I: 
Labbe et al.l 120061 : iGrazian et al.l 120061 : Ide Santis et al.l 
2007ft . Here, we introduce basic concepts of TFIT to 
highlight the main difference between this work and 
previous ones in the literature. We further demon- 
strate the effectiveness of TFIT for faint high-redshift 
galaxies by directly comparing the TFIT results with 
the conventional methods, namely, aperture photometry 

The basic premise of TFIT is that the best flux of a 
source can be estimated by simultaneous fitting of all the 
adjacent sources rather than just performing photome- 
try within a fixed aperture. This strategy is particularly 
effective for crowded fields where sources are often sub- 
stantially blended with one another and PSFs are large 
enough to not just affect their immediate neighbors but 
possibly those further out from any given source. Even 
though deep fields such as GOODS are by design mostly 
devoid of bright sources in the optical bands, the same 
field is crowded in the deep IRAC data. 

In order to achieve simultaneous fitting of multiple 
sources, TFIT constructs a realistic model of how an 
object would appear in a low-resolution data based on 
the observed morphology of the object in a deep high- 
resolution image. The cutout image of each object is 
created from the high-resolution data, then convolved 
with the kernel designed to reproduce the PSF of the 



low resolution image and block-averaged to the desired 
pixel scale. The image created with this procedure is 
called a template, which provides the best approximation 
of how the same object appears in the low-resolution im- 
age. This assumes no significant dependence of galaxy 
morphology on wavelength - i.e., no morphological k- 
correction. While such an assumption may not be valid 
in some cases, this can be in part overcome by using the 
high-resolution image observed in the wavelength range 
closest to the low-resolution one. For the GOODS TFIT 
photometry, we use the Zsso band to create templates 
for IRAC photometry and the S435 band for the U-b&nd 
photometry 

Once the templates of all galaxies are created, TFIT 
performs chi-square minimizations on subsections of the 
data, to find the best- fit fluxes for each object. Fur- 
thermore, TFIT at this stage corrects for any second- 
order astrometric misregistration which may be present 
for data sets taken with different instruments. Such mis- 
registration can bias the measured flux and colors not 
only for TFIT but also for aperture photometry as well. 

We carried out simulations to compare the perfor- 
mance of TFIT in direct comparis ons with aperture 
photo metry returned by SExtractor (|Bertin fc Arnoutsl 
1996). Artificial galaxies (with known input fluxes) are 
inserted into the real data (ACS and IRAC) and fluxes 
are recovered by both aperture photometry and TFIT. 
As for aperture photometry, we adopt MAG_AUTO for 
both ACS and IRAC as the best-estimate. Figure Q] 
shows the z — [3.6/xm] colors measured from SExtractor 
(left) and TFIT (right) in comparison with the true 
input colors of the same sources. It is evident that 
the colors measured by TFIT are much closer to the 
true values than those by aperture photometry, and 
also that photometric scatter is considerably smaller. 
Furthermore, catastrophic outliers identified in aperture 
photometry (shown in red in both panels) are minimized 
in TFIT photometry as they are located much closer 
to the one-to-one line (green) in TFIT colors. This 
is not surprising because such errors typically result 
from significant contamination due to crowding. Such 
catastrophic failures are expected to increase in deeper 
images where the source density is higher. We also note 
that such cases are expected to be frequent for faint 
high-redshift sources whose flux can be significantly 
affected by usually-brighter neighbors when the conven- 
tional method is adopted. 

Throughout this work, we utilize the GOODS Grand 
Unified TFIT catalo£3 to measure the spectral proper- 
ties of LBGs. The catalog includes the SExtractor pho- 
tometry of the HST /ACS band (-B435, ^06, 1775, z S50 ) 

11 The catalog is soon to be made a vailable at 
http : //www. stsci . edu/science/goods/DataProducts/ 
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Fig. 1. — Optimal multiwavclcngth photometry by TFIT. Extended sources with a range of colors are simulated to compare the perfor- 
mance of conventional aperture photometry and TFIT. (Left) The z — chl (ACS Z850-band - Spitzer [3.6/xm]) colors measured by SExtractor 
are compared with the input colors for 600 simulated galaxies. Large scatter is expected due to source blending and confusion. (Bight) 
The colors of the same sources measured by TFIT. The sources with the "aperture color" bias >1 mag are shown in red on both panels. 
Note that most of the same sources are successfully recovered in TFIT-derived colors. 



and TFIT photometry for the remainder of our pho- 
tometric bands (UJHKs [3.6^m][4.5^m][5.8^m]). The 
ACS zg5o-band (F850LP) was used as the detection band 
for the source extraction and to define isophotal area 
for other ACS bands, and also as high-resolution tem- 
plates for the TFIT photometry. This ensures that the 
measured flux either by SExtractor or TFIT will re- 
turn accurate colors within the same isophotal apertures, 
and that aperture correction is straightforward to han- 
dle. Galaxy colors are computed using the MAGJSO 
for the ACS bands and TFIT fluxes for the others, while 
MAG_AUTO in the zgso-band is taken as the total mag- 
nitude. The limiting magnitude for each band is esti- 
mated from the distribution of TFIT photometric errors 
(la) for all galaxies in the catalog, which typically ex- 
hibits a well-defined peak. The peak value of the error 
distribution represents the typical uncertainty given the 
sensitivity of the data. Whenever TFIT fluxes are con- 
sistent with zero within errors, we take the la errors as 
upper flux limits. 

Finally, we have compared the GOODS TFIT cata- 
log with t he public GOODS-M USIC catalog (v2.0) dis- 
cussed in San tini et al.l ((2009) , which uses a similar 
algorithm (CONPHO T) to deal with source blending 
(|de Santis et al.l 120071 ) . The photometry from the two 
catalogs are found to be consistent with each other with 
the median offset less than 0.1 mag down to 26.0 mag 
(for IRAC [3.6/im]). The scatter, however, is quite large 
(« 0.27 mag) towards the faintest bin. While the origin 
of the relatively large scatter is beyond the scope of this 
work as it is difficult to quantify without detailed anal- 
yses coupled with simulations, no systematic bias found 
between the two catalogs implies that our main results 
will remain robust against the particular choice of algo- 
rithms adopted for optical photometry discussed above. 

3. MEASURING STELLAR MASSES OF HIGH-REDSHIFT 
GALAXIES 



3.1. SED Fitting 

Stellar population synthesis models are used to derive 
stellar masses of galaxies in our sample. We use the up- 
d ated version of the stellar population synthesis models 
of lBruzual fc Charlotl (|2003l known as the BC03 model) , 
which we refer to as the CB07 code hereafter. We use 
the CB07 models to compute the broadband color evo- 
lution of galaxies wi th differen t star formation histories 
and ages. We use a IChabrierl (|2003l ) initial mass func- 
tion (IMF) with a lower and upper mass cutoff at O.IMq 
and 100M Q , respectively, while the metallicity is fixed 
to solar. While the CB07 model includes enh anced con- 
tribu tion from the TP-AGB stars (see, e.g., iMarastor] 
2005), it returns stellar masses very similar to the BC03 
because the light from UV-selected galaxies is dominated 
by younger OB stars even at the rest-frame optical wave- 
lengths sampled by the IRAC channels at z > 3. The 
TP-AGB stars would contribute more at longer wave- 
lengths (rest-frame near-infrared) and for older stellar 
popu lations (> 1 Gyr, also see Figure 5 in Sta rk et al.1 
[2009h . 

The model galaxy SEDs are redshifted in the range of 
z = 0.2 — 8 with Az = 0.1. Internal dus t extinction is 
compu ted according to the prescription of iCalzetti et al.l 
(2000) law in the range of E(B - V) = 0.0 - 0.95 
with AE(B — V) = 0.025. We also attenuated the 
resulting SEDs from the neutral hydro gen abs o rption 
in the intergalactic medium using the iMadaul (|1995» 
prescription. The age of our model SEDs extends 
from 5 Myr to the age of the universe at the given 
rcdshift. The star formation history is parametrized as 
an exponentially decreasing star formation rate (SFR), 
where r represents the e-folding decay time. We use 
r = 0.1,0.2,0.3,0.4,0.6,0.8, and 1.0 Gyr, and constant 
star formation (CSF; r — > oo). We precompute a large 
number of model SEDs spanning the parameter space 
described above, then integrate them through the filter 
response functions of all the observed bands (HST / ACS: 
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-B435^606i7752850, Spitzer/IRAC: [3.6/mi], [4.5/xm], 
[5.8/Ltm], VLT /ISAAC: JHK S , CFHT/WIRCAM: JK S , 
KPNO/MOSAIC L/-band, and VIMOS/fJ-band) to 
compute the galaxy colors for each point in our model 
grid. 

Using the precomputed model grid, we find the best- 
fit parameters for each galaxy in our sample. To evalu- 
ate the goodness-of-the-fit, we first determine the overall 
normalization to match the observed photometry, then 
compute the chi-square using all relevant data points. 
Because model SEDs are normalized to have total mass 
of IMq, the normalization factor directly determines the 
stellar mass. For each galaxy, we also determine the ab- 
solute UV magnitude at 1700 A, Mn 00: by interpolat- 
ing the best-fit model spectra. The evaluation of the 
chi-square value is done in flux units rather than mag- 
nitudes, which helps properly handling the data points 
with non-detection. Whenever the photometry in a given 
passband is consistent with zero within the la error, we 
assign as a flux upper limit a value that is representative 
of the sensitivity limit for the band. This treatment is 
often useful to avoid a large stellar mass to be assigned as 
a solution when TFIT returns very large photometric un- 
certainties. Such cases can occur when a source is hope- 
lessly blended with much brighter neighbors. However, 
in most cases, using TFIT-returned flux and flux errors 
(instead of flux upper limits) for the fitting should yield 
similar results. Data points with upper limits are in- 
cluded in the chi-square evaluation only when the model 
flux is higher than the upper limit. 

The flux upper limits in these bands are listed in Ta- 
bl e [H which a r e simi lar but not identical to those given 
in IStark et al.l (|2009f) for the same data. The main rea- 
sons for this slight difference are two-fold. First, we used 
the source-weighted RMS map rather than the exposure 
map, which helps make more realistic estimate of the 
noise near bright sources. Second, TFIT includes the 
flux uncertainties arising from source blending, which, in 
particular, for faint sources like most high-redshift galax- 
ies, is an important factor to account for as the majority 
of them are blended with brighter neighbors. 

Instead of using the conventional approach of finding 
a single value as the best-fit stellar mass, we determine 
the full probability distribution function for each galaxy. 
In practice, this means that we retain all models that 
are within Ax% = 1 from the minimum \ 2 model as 
acceptable fits, and thus the stellar mass of the galaxy 
is described by a histogram of stellar masses spanned 
by all acceptable models. While this approach return 
results that are fully consistent with the conventional 
method when stellar mass of a given galaxy is in ques- 
tion, it has several advantages from a statistical stand- 
point. For example, there can be multiple models with 
X 2 values very close to one another but over a range of 
stellar mass. These situations become very real when 
one is dealing with galaxies at the low-mass end. Pick- 
ing one model over the others may bias the outcome in an 
unpredictable fashion unless the systematics that affect 
the fitting are well understood. Choosing to construct 
the number counts using the full probability distribution 
makes it easier to correct the "observed" number counts 
for the measurement biases and systematics, which can 
be quantified using extensive simulations of galaxies. In 



§4, we present our simulation procedures in detail, and 
further discussions on how best to measure the galaxy 
counts in stellar mass bins given the photometric uncer- 
tainties. 

The redshift distributions and the accuracy of pho- 
tometric redshifts for our samples are w ell understood 
(jVanzella et al.|[200l iDahlen et al.ll2(H0l ). Thus, we rely 
on photometric redshift estimation to derive physical pa- 
rameters of each galaxy. During the SED fitting, red- 
shift is varied within the best-fit photometric redshift by 
±Az = 0.125, which accounts for a typical uncertainty 
of the photometric redshift estimation (see Figure 9 in 
IDahlen et al.l [20101) . except when the accurate spectro- 
scopic redshift is known. Varying the redshifts is im- 
portant in obtaining realistic errors in the derived SED 
parameters for these galaxies as the majority (95% of the 
full i?435-band dropout sample) does not have a secure 
spectroscopic redshift. Typically, the errors in stellar 
mass when the redshift is varied are 20-25% larger than 
those when the redshift is fixed to the best-fit photomet- 
ric redshift. 

3.2. The Observed Number Counts of Galaxies 

Following the procedure described above, we measured 
the stellar masses of the -B435- and Vgo6-band dropouts 
and constructed the number counts. In order to assess 
robustness of our results, we investigate how the final 
number counts depend on the assumed star formation 
history. In Figure [2] we compare the number counts in 
the stellar mass bins for the i?435-band dropouts when 
different star formation histories are assumed. While we 
only include constant star formation history (CSF) and 
exponentially declining r models in our fitting procedure, 
comparing the results using a short r model (100 Myr, 
in this case) and CSF model gives an idea of how the 
results would differ for scenarios in which star formation 
is assumed to take place on short, "bursty" timescales, 
versus longer, more continuous modes. 

We find that the observed number of galaxies in a given 
mass bin is slightly higher when a larger r value is as- 
sumed as shown in the left panel of Figure [21 This is be- 
cause models with a larger r typically result in a higher 
stellar-mass-to-UV-light ratio (mass-to-light ratio, here- 
after), and thus predict slightly larger stellar mass for 
given rest-frame optical fluxes. However, different mass- 
to-light ratios seem to result in a solid shift rather than 
change in the overall shape of the number counts. They 
agree well with one another when a small offset is applied 
to correct for different mass-to-light ratios. In the right 
panel of Figure [H we show the same number counts with 
a small offset A log[M star /M ] = 0, 0.13, 0.03, 0.0 for the 
free-r model, r=100, 300, and CSF model, respectively. 
This shows that the overall shape of the number counts 
can be measured robustly regardless the assumed star 
formation history. One possibility is that a typical r 
value depends systematically on galaxy UV luminosity or 
stellar mass. For example, if low-mass galaxies have, on 
average, shorter SF timescales, such a dependency would 
manifest as a slightly steeper slope on the low-mass end. 
While it is difficult to model such an effect without more 
observational or theoretical input, Figure suggests that 
such an effect should have a relatively small impact on 
the galaxy counts. Throughout this work, we use the 
number counts measured with "free-r" models for our 
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Fig. 2. — Left: The observed number counts of LBGs at z ~ 4 in stellar mass bins when different star formation histories are assumed. 
Larger r values generally return higher stellar-mass-to-UV-light-ratios, and thus slightly larger masses for given rest-frame optical fluxes. 
Right: When a small offset is applied to the masses to correct for this effect, the number counts using different SFHs agree with one 
another very well. This means that even if mass measurements are biased towards lower or higher-masses because of our assumption in 
star formation histories, the overall shape of the number counts is likely robust. 



analyses unless stated otherwise. 

There are several factors that may contribute, to a 
varying degree, to the estimation of the number counts. 
First, a hidden population of old stars formed from a pre- 
vious generation of star formation may contribute to the 
stellar mass. The upper limit on this contribution can be 
set by assuming a population formed from an instanta- 
neous burst at a very high redshift, which maximizes the 
mass-to-light ratio. Based on the multi-wavelength data 
set (out to ii'-band, which corresponds to the rest-frame 
U-band) for z ~ 3 g alaxies selected from the HDF-N, 
IPapovich et al.l (|2001|) estimated that galaxies, on aver- 
age, can be hiding stellar mass up to a factor of 3. While 
this can certainly be the case for some galaxies in our 
sample, the relative contribution from an old popula- 
tion is likely smaller because the universe is younger at 
z ~ 4 and ~ 5, and thus, the galaxies have had less 
time to have formed a significant stellar population un- 
related to the current star formation. Furthermore, the 
Spitzer data samples out to, at least, the rest-frame z 
(J)-band, «9000 (7500) A, at z = 4 (5) as the majority 
of the galaxies with stellar mass M sta r > 10 8 ' 7 M Q are 
detected in the [3.6/Ltm] and [4.5/xm] bands. For galax- 
ies which are not formally detected in the Spitzer bands, 
it is difficult to assess how much mass may be hidden 
in old stellar populations. However, the average SED of 
Spit zer-iaint galaxies c onstructed by stacking (see Figure 
8 of iStark et al.l [20091 ) suggests that it is unlikely that 
such a contribution is more substantial than their more 
luminous counterparts. Furthermore, a large population 
of maximally old galaxies would be inconsistent with the 
rapid redshift evolution of the UV star formation rate 
density and UV luminosity function observed at z > 4 
(jBouwens et al.|[2007l I2008T) . 

Another inte r esting possibility concerns a recent result 
by iShim et all (|2011l ) , who reported a strong Ha line 
emission in a large fraction (>70%) of galaxies detected 
in the Spitzer [3.6/mi] and [4.5/im] with high signal-to- 
noise ratios. While the study is based on a relatively 
small sample (70 galaxies) of spectroscopically confirmed 
sources at 3.8 < z < 5.0, they convincingly showed that 
the flux excess in the Spitzer [3.6/im] band is indeed due 



to strong Ha emission from star formation and not from 
hidden AGN activities. From the viewpoint of measuring 
stellar mass using broadband photometry, the presence 
of strong emission line can skew the mass estimation be- 
cause emission lines are not included in our SED model- 
ing. Furthermore, the observed equivalent widths of the 
Ha emission, they claimed, suggest that it requires ex- 
tremely young ages, top-heavy IMF, or sub-solar metal- 
licity (or a combination of these). In practice, IRAC 
[3.6/im] is one of the main drivers in the SED fitting 
procedure that determines the stellar mass as it samples 
the rest-frame /-band. We visually inspected the ob- 
served SEDs of the i?435-band dropout sample to check 
whether such an excess is common also in the photo- 
metric sample. While the effect of line emission in the 
photometric data is fairly subtle to identify visually, and 
may be particularly hard to see for the fainter galaxies 
without spectroscopic redshifts, we determine that the 
number of sources with a noticeably large [3.6/im] ex- 
cess is indeed small, and thus unlikely result in change 
in the overall shape of the galaxy number counts. Be- 
cause the strength of Ha emission line should be lower 
for UV-faint (l ower SFR) galaxie s (a typical zgso-band 
magnitude for Shi m et al.l (|2011[) sample is 25.0, while 
the majority (> 95%) of the galaxies in our sample are 
fainter than 25.0 mag), the lack of Ha-excess sources 
seems to be in line with a possibility that the line emis- 
sion is indeed due to star formation. Regardles s of th e 
nature of the flux excess reported by IShim et al.l (2011), 
we conclude that it is unlikely that it will change our 
results significantly. 

3.3. The Observed Muv-M sta . r Scaling Law 

Having measured both UV luminosity and stellar mass 
of the dropout samples, we present how these galax- 
ies populate the plane of the observed luminosity and 
stellar mass in the main panels of Figure [3] The fig- 
ure is color-coded by the total number of galaxies in 
each given cell, and thus illustrates the overall trend of 
the relative location of galaxies on the plane as well as 
the number of galaxies at different mass and luminosity 
bins. The average stellar mass at a fixed UV luminos- 
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ity is also computed from our measurements as shown 
as yellow filled circles. Top and right panels indicate 
the sum of all columns and rows (of the main panels), 
i.e., galaxy counts as a function of UV luminosity and 
stellar mass, respectively. We also note that a large 
concentration of UV-faint galaxies at very low masses 
(10 8 - 10 8 - 5 M Q ) is partly an artifact produced by a 
combination of flux upper limits and our adopted fit- 
ting method. When a source is not formally detected in 
the Spitzer IRAC channels, the minimum x 2 is achieved 
when the model flux is below all the flux upper limits set 
by the data. As a result, the best- fit solution is returned 
with a range of low stellar masses, which is primarily 
determined by a range of stellar-mass-to-UV-light ratios 
allowed within our model grid (including the models with 
the lowest mass-to-light ratios as the data does not pro- 
vide any constraints), and thus does not reflect a true 
value. To avoid any systematic bias in our results by 
this effect, we confine our analyses only to galaxies with 
stellar mass M star > 1O 8 - 7 M for the £?435-band dropouts 
and M star > 10 9 M© for the Vej06-band dropouts. These 
limits, the minimum stellar mass considered in our anal- 
yses, are also indicated in Figure [3] on right panels as 
"Minimum Mass". Our simulations (S|4|) indicate that 
the current data allow statistically robust measurements 
down to this limit. A majority of simulated galaxies in 
this mass range are formally undetected in the IRAC 
bands, or have large photometric errors. 

Our measure s are l argely consistent with those made 
by iStark et "all (|2009f ) marked by dashed lines, but are 
generally slightly lower (by ~ 0.1 dex) than their esti- 
mate. The discrepancy may be due to the fact that they 
adopted a single value for each stellar mass rather than 
taking into account the full distribution, and as a result, 
the median is slightly skewed towards higher masses. As 
we will show in £j4j the uncertainty in the mass determi- 
nation at 10 9 Afo is somewhat large, and thus requires 
statistical approach to construct the distribution. Fur- 
thermore, the discrepancy may reflect the difference in 
the two samples, that ours includes a much larger num- 
ber of galaxies («90% of the photometric sample in com- 
parison to their 35%) because TFIT photometry does not 
require that sources be well isolated. 

In Figure [H we show that there is an overall correla- 
tion between stellar mass and UV luminosity down to 
Mstar ~ 1O 8 ' 7 M0 in that more UV-luminous galaxies are 
on average more massive. However, it is equally evident 
that there is a substantial scatter to this correlation. In 
particular, there exist galaxies that are quite massive but 
faint in the UV. Roughly 36% and 20% of the galax- 
ies in the -B435- and V(306-band dropout sample that are 
more massive than Af sta r = 10 9 5 M Q are fainter than 
A/1700 = —20. The implication is that there is a non- 
negligible population of galaxies whose SFRs were con- 
siderably higher in the past, but are currently fading in 
the UV. It is possible that some of them may be faint sim- 
ply due to more dust reddening, but as we will show later 
( H6.ip . dust alone is not enough to explain the observed 
UV luminosity distribution of these low-mass galaxies. 
Interestingly, the opposite is not true in that the region 
of the low stellar mass and high UV luminosity is largely 
devoid of galaxies. This suggests that the majority of 
the UV-luminous galaxies have been forming stars for 
quite some time (at least a few hundred million years). 



This interpretation is consistent with a recent study of 
the most UV-lu minous galaxies at the same redshift by 
iLee et all (|2011l also shown in Figure [3]) who concluded 
tha t they undergo sm ooth star formation histories (also 
sec, IShim et alJl20ll . 

4. UNDERSTANDING MEASUREMENT 
UNCERTAINTIES/BIASES 

Spectral energy distribution (SED) fitting is a proce- 
dure that attempts to fit a number of galaxy parameters 
based on multiple data points taken at different depths 
and resolutions. Thus, its outcome may depend strongly 
on the range of galaxy parameters (initial mass function, 
star formation histories, age) and the quality of the pho- 
tometric data. Therefore it is worth paying close atten- 
tion to how well each galaxy parameter can be measured 
given the set of observations. The most relavant issues 
to our analyses are: How robust are the derived galaxy 
properties (e.g., extinction, stellar mass) when the SFH 
and redshift are unknown? What are their uncertainties 
and how do they scale with the quality of photometry? 
How well can we measure the galaxy light within the 
same effective aperture from the multi-resolution data 
set, using a photometric technique of choice (e.g., SEx- 
tractor, TFIT)? Creating a mock data set with realistic 
properties (noise, seeing, galaxy colors, etc.) is key to 
properly accounting for all factors that may contribute 
to the measurement errors and systematics. To address 
these questions, we carried out a large set of simulations. 
The procedure is designed to simulate galaxies as realis- 
tically as possible in every step of the way by mimicking 
the noise properties, size distribution, photometric ex- 
traction, and derivation of physical parameters via SED 
fitting, of the real high-redshift galaxies. 

4.1. Galaxy Simulations on Mixed Resolution Dataset 

The simulation procedure can be summarized as fol- 
lows: 1) Create galaxy model SEDs spanning a wide 
range of physical parameters (age, reddening, SFHs, red- 
shift, and stellar mass), 2) Compute input photome- 
try of these SEDs in the observed passbands, 3) Insert 
mock galaxies into the images in the different photomet- 
ric bands, 4) Measure photometry using TFIT software, 
and 5) Carry out SED fitting to "measure" output physi- 
cal parameters. By design, the output of our simulations 
is identical to that of the real galaxies except that input 
parameters are known for simulations. By comparing the 
input and output parameters, we can quantitatively as- 
sess the robustness of the measurements at various levels 
and use that information to correct for the systematics. 
As for the range of physical parameters used to create 
mock galaxies, we use the same ranges that we explore in 
the SED fitting procedure for real galaxies. By doing so, 
we ensure that we can quantify all the measurement er- 
rors and biases within the assumptions that went into de- 
riving the physical parameters of real galaxies. Any non- 
standard physical characteristics of real galaxies, such as 
exotic IMFs or sub- or super-solar metallicities, that may 
have contributed to our measurements, however, cannot 
be accounted for in our simulations, and thus is beyond 
the scope of this work. The descriptions of the steps #1 
and #2 are given in i j3.ll The details of inserting mock 
galaxi es in the real data (#3) are found iFcrguso n et al.l 
(|200l . 
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Fig. 3. — The number counts of LBGs in UV luminosity and stellar mass bins at z ~ 4 (left) and z ~ 5 (right). We show in the main 
panels, the location of all galaxies in our sample. Each cell is color-coded with the observed number of galaxies. The average stellar mass 
at a fixed UV magnitude is also shown as yellow filled symbols. Our results indicate an overall correlation betwe en UV luminosity and 
stellar mass albeit with a significant scatter. Our estimates are in good agreement with those of Stark ct al. (2009) as marked by dashed 
line and fi l led sy mbols, but extend further down to lower masses using TFIT photometry. For z ~ 4, we also show similar measurements by 
Lee ct al. (2011) for very UV-luminous (Mi 700 < — 21) galaxies. In upper panels and right panels, we show the galaxy number counts in UV 
magnitudes and stellar mass bins, respectively. The minimum mass we include in our stellar mass function analyses, which we determine 
based on our s i mulat ions (i|4t. is also indicated on right panels. Our number counts in UV magnitudes are in excellent agreement with 
Bouwcns ct al. (2007), confirming that the UV luminosity function has a steep faint-end slope. A large concentration of UV- faint galaxies 
in low masses is due largely to non-detection in IRAC channels and thus represents upper mass limits in some galaxies. This suggests that 
at Mi 700 > — 19 the median stellar mass for galaxies falls more steeply than that for more luminous ones. 



The main improvement of our simulations over the 
previous version described in iFerguson et al.l (|2004f ) is 
that we extend the technique to create multi-resolution 
data sets. This revision is straightforward as all of the 
data have pixel scales that are integer-multiples of one 
another (O'.'03/pix for the H ST /ACS data, O'.'5/pix or 
O'.'30/pix for the ground-based data, and O'.'60/pix for 
the Spitzer /IRAC data). All images are initially created 
in the highest resolution, convolved with the appropri- 
ate PSFs, then block-summed to the desired pixel scale 
before being added to the real science images (step #3). 
Only 30 galaxies are inserted at a time over a 4 arcmin 2 
area (4096 x 4096 in the pixel resolution of the ACS 
bands) to avoid unrealistic crowding. They are placed at 
random locations within the images with random orien- 
tation angle. The light profiles of galaxies are chosen to 
be either of De Vaucouleurs profile or of exponential disk. 
The size distribution is set to follo w a lognormal distribu - 
tion as observed at these redshifts ((Ferguson et al.ll2"004f ). 

For each run, once the mock images are created, we 
repeat the identical procedures for measuring their pho- 
tometry (step #4) and SED parameters (step #5) as for 
real galaxies. This means that we run SExtractor on the 
ACS images with the same setup for source detection and 
photometry, then use the ACS zsso-band as templates to 
run TFIT on other passbands. The J7-band (below the 
Lyman limit for z > 3 galaxies) was not included in the 
simulations, as mock galaxies do not suffer from low- 
redshift interlopers, for which the L/-band is mainly used 
to check non-detection for real galaxies. 



4.2. How Well Can We Measure Galaxy Parameters': 



Using the procedures described above, we have simu- 
lated about 200,000 galaxies, of which «35,000 galaxies 
are selected as £?435-band dropouts and «26,000 as V606- 
band dropouts using the color selection criteria given in 
Equation 1 and 2. These galaxies span a wide range 
of UV luminosities, star formation histories, age, extinc- 
tion, and stellar mass. In Figure HI we show the input 
and recovered value of population age, dust reddening 
(parametrized as E(B — V)), and stellar mass, for simu- 
lated galaxies. Stellar mass is the most robust quantity 
that we can recover based on the available broadband 
photometry at the current depth, while internal dust 
reddening can be recovered to a somewhat lesser degree 
(bottom right). As for population ages, our SED fitting 
fares very poorly with a large scatter (top right) unless 
the age of input galaxies is much older (> 500 Myr) than 
the values derived by SED fitting for the majority of the 
observed Lyman break galaxies. The poor recovery of 
population age is likely attributed to the current (shal- 
low) depth of the near-infrared data from ground, which 
samples near the Balmer/4000A break (H and K$ band 
for the i?435-band dropouts at z ~ 4). These resu lts are 
in qualitative agreement with iLee et all (|2009bO . who, 
based on the SED fitting of mock galaxies created by a 
semianalytic model, reported that stellar mass can be re- 
covered most reliably (within 0.1 dex of the true values) 
even when the assumed star formation history is clearly 
different from the intrinsic one. 

In Figure [5j we show typical SEDs of real (left) and 
simulated (right) galaxies at different stellar masses. For 
the simulated galaxies, we also indicate the input pho- 
tometry as open circles, illustrating the increasing de- 
gree of photometric scatter and errors towards the low- 



10 



o.o 




1.0 -0.5 0.0 
log[Age ;n (Gyr)] 



§ 10 



en 
o 




0.2 0.3 0.4 0.5 
E(B-V) ;n 



Fig. 4. — The probability distribution function for the stellar mass (lower left), internal reddening (lower right), and age parameters 
(top right) recovered from our SED fitting procedure as a function of input parameters. The color bar (top) indicates the probability value 
for each (input, output) cell. The discreteness of input ages reflects non-equal steps of ages adopted for the simulated galaxies. These 
results are based on 34,501 f?435-band dropouts planted in the GOODS-S and GOODS-N field according to the simulations procedure we 
described in text (33). Based on these results, we determined that stellar mass is the most robust quantity that can be recovered, at least, 
based on the broad-band photometry via SED fitting. The internal reddening can be recovered to a somewhat lesser degree, indicated by 
a larger error. The population age is the least robust parameter, which is likely driven by the current sensitivity limit of the near-infrared 
data sampling the Balmer/4000A region (H and Kg band for the £?435-band dropouts). 



mass end. On right panels, we show the probability dis- 
tribution of the recovered stellar mass at various input 
mass values. As expected, it is evident that the accuracy 
of stellar mass estimation decreases rapidly towards low 
masses as a result of larger photometric errors in rele- 
vant passbands (in this case, mainly driven by the IRAC 
[3.6/im] and [4.5/im]). Increasingly large spread of recov- 
ered mass with respect to the true value further supports 
the statistical approach we take in constructing galaxy 
counts in mass bins. In any stellar mass function that 
rises relatively steeply towards the low-mass end, such 
scatter will result in artificial steepening of the observed 
slope as more number of low-mass galaxies will scatter 
into higher-mass bins than the opposite. 

In our analyses throughout this paper, we use the 
mass probability distribution (right panels of Figure [5]) 
for two purposes; first, to assess and thereby determine 
how much we can push the limit of the available data 
without introducing systematic biases into our measure- 



ments, second, to statistically correct the observed num- 
ber count for the "mass spread" due to photometric er- 
rors. As for the former, we conclude that stellar mass can 
be measured reliably down to log[M stal /Mo} = 8.7 and 
9 at z ~ 4 and 5, respectively. As evident in right panels 
of Figure [3 the probability distribution becomes skewed 
towards the masses below the limit. This is mainly due 
to the fact that sources become very faint (generally, in 
all bands). This leads to two independent consequences 
that affect the SED fitting procedure. First, the SEx- 
tractor tends to underestimate the isophotal area from 
the detection image, which is later used to define the 
"effective aperture" in other bands to measure flux using 
TFIT. More quantitative discussions on this issue will be 
discussed in an upcoming paper (A. Galamctz ct al., in 
prep). Second, because the same sources are (usually) 
intrinsically fainter in the Spitzer bands, and as a result, 
the SED fitting procedure is often driven by the flux up- 
per limits. Because there is no meaningful constraint 
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on the rest-frame optical wavelengths, the fitting returns 
a wide distribution of stellar masses as acceptable (see 
discussions in H3.2\i . 

5. THE UV LUMINOSITY FUNCTION AND STELLAR 
MASS FUNCTION OF LBGS 

5.1. The UV Luminosity Function 

We measure the UV luminosity function of LBGs in our 
sa mples using a method ology similar to that presented 
bv iBouwens et alJ (|2007h . The procedure consists of 1) 
constructing the number counts of the real galaxies in 
apparent magnitude bins, 2) quantifying the photomet- 
ric and selection completeness of the observed galaxies 
by using simulations (Q, and 3) determining the range 
of UVLF parameters that reproduces the observed num- 
ber counts. In quantifying the completeness, we treat 
internal dust reddening as an explicit variable. 

We define the selection efficiency, ps{m' , z,e), as the 
probability that a galaxy at redshift z of the (true) 
apparent magnitude to' with the reddening parameter, 
E(B — V) = e, to be detected and selected as a dropout 
in our selection. In addition to the selection efficiency, 
we further quantify the photometric error that affects the 
number counts of the dropouts as the probability that a 
galaxy (m',z,e) to be observed at the apparent magni- 
tude to, hereafter referred to as pm(w|to', z, e). To re- 
late the intrinsic UV luminosity function to the observed 
number density, we define the effective volume as: 

f C d\^ 

V m , k = / W(M(m',z,e) - M k )V{m\m',z,e)—dm'dz 

~ zZYl W(M(m', z, e) - M k )V(m\m', z, e )^( z ) 

m' z 

xAra'Az (1) 

where W is a window function, and V(m\m' , z, e) = 
Ps(m, z, e) x pM(m\m', z, e), dV/dz is the cosmic volume 
covered by the survey area, to is the observed apparent 
magnitude, to' is the true apparent magnitude, and M k is 
the absolute UV magnitude at 1700 A. The effective vol- 
ume averaged over the intrinsic distribution of reddening 
e is: 

,,r v _ ^ n V m ,kT>(e n , M k )Ae 
E n T>(e n ,M k )Ae 

where T> characterizes the distribution of reddening 
which can be luminosity-dependent (M k ). The observed 
number count for a given luminosity function can be com- 
puted as: 

N m = Y,MM k ){V m ,k)e (3) 

k 

where 0i(Mfc) is the UV luminosity function evaluated 
at the absolute magnitude bin Afi7oo = M k , where the 
UVLF is parameterized as a Schechter function: 

^(M fe ) = (0.41nl0)<f io°- 4 ( m uv-a^)(i+«) 

xexp[-10 - 4 ( M «v-^))] ( 4 ) 

As for the reddening distribution £>, we assumed a nor- 
mal distribution with a mean of e = 0.14 if M k < 
-21.06 (£*), e = if M k > -18.56 (0.1X*), and linearly 



declining in between. The la scatter of the distribu- 
tion remains o~ e = 0.14 at all luminosities . Our prescrip- 
tion is essentially identical to that used bv lBouwens et all 

In Figure [6l we show the ob served number cou nts to- 
gether with those published in iStark et al.l (|2009T > based 
on the GOODS vl.O catalog. Both measurements agree 
well except in a few faintest magnitude bins, where our 
deeper catalog should be more complete. Thick dashed 
lines indicate the expected number counts computed ac- 
cording to E quation [3] when the be st-fit UVLF parame- 
ters given bv lBouwens et al.l (|2007l ) are assumed. While 
we are unable to independently determine the UVLF 
parameters based on the GOODS data alone, the good 
agreement between the two suggests that the faint-end 
slope is steep in the range of —(1.7 — 1.8). To illus- 
trate the fact that our data is indeed more consistent 
with a steep (< —1.7), we also show in Figure [5] as a set 
of dotted lines the model number counts with the slope 
— (1.6,1.4,1.2) when the characteristic luminosity M{j v 
is fixed. All functions are normalized at z^q — 24.97 and 
25.71 at z ~ 4 and 5, roughly corresponding to the re- 
spective characteristic luminosity at the median redshift 
of each sample (shown as filled arrows) . 

In comparison with the results by IBouwens et al.1 
pool , we find a slight discrepancy at the bright end, 
where their es timate is about a factor of 2 higher than 
what we and IStark et all (|2009D found. This is likely 
due to the fact that we discarded some number of bright 
sources after inspecting their SEDs (as well as MIPS de- 
tection, see H2.3P . Many of bright sources have UV colors 
consistent with high-redshift LBGs (and thus included in 
their sample), but longer wavelength data suggests that 
they are unlike ly to be at high r edshif t. Clearly, both our 
estimates and IBouwens et al.l (|2007| ) suffer from small- 
number statistics at the bright end, and the discrepancy 
is driven by mere 10-20 sources and their photometric 
redshifts, more complete spectroscopy of these sources 
and larger-area surveys will be needed to better quantify 
the UVLF at the bright end. 

5.2. The Stellar Mass Function of SFGs 

Using robust measurements (and error estimates) of 
SED properties and UV luminosities for individual galax- 
ies, we next determine the stellar mass function of LBGs 
at z ~ 4 and 5. We caution readers that these mea- 
surements do not in any way represent the total stellar 
mass function of all galaxies, but rather provide informa- 
tion on how much of the cosmic stellar mass budget is 
contributed by actively star-forming and relatively UV- 
bright galaxies. Furthermore, as we will see later, by 
comparing the shape of the SMF with that of the UV LF 
measured for the same galaxies, one can obtain crucial 
clues to how (UV-visible) star-formation has proceeded 
with cosmic time. 

We compute the number density of LBGs at each stel- 
lar mass bin as follows. First, the number density of 
galaxies at a fixed apparent magnitude rnj and stellar 
mass M* jm is obtained by correcting the observed num- 
ber N b s (mj,M* im ) for the selection incompleteness by 
dividing it by the effective volume V e s(mj). Using Equa- 
tion[T]and[2l the total effective volume at a fixed apparent 
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Fig. 5. — We show the spectral energy distributions of typical galaxies at different stellar mass bins for both real (left panels) and 
simulated galaxies (middle panels). Filled circles represent flux in each photometric band with TFIT errors, while solid lines indicate the 
best-fit population models from which the stellar mass is determined. All stellar masses are given in units of solar mass. In the histograms 
on right panel we show the distribution of recovered stellar mass with respect to the input mass. The typical accuracy in the determination 
of stellar mass can be assessed based on the mean and dispersion of the output mass distribution. We conservatively set the mass limit to 
log [M star /MQ] = 8.7 and 9.0 at z ~ 4 and 5, respectively. Below these limits, output masses are increasingly uncertain as evidenced by 
an offset with large scatter (top right panel). 



magnitude rrij is 



(t> M (mj,M.. 



_ N ohs (m,j,M* tm ) 



V e s(mj) 



(5) 
(6) 



The total number density of LBGs at a stellar mass bin 
is obtained by adding the contributions from all the ap- 
parent magnitude bins: 



(7) 



To estimate realistic errors at each mass bin, we ac- 
count for both random Poisson errors (based on the ob- 
served number at each bin) as well as sample variance 



(based on the difference between the observed number 
at two fields). The Poisson error at each bin is com- 
puted as A^j fe f(m i j,M* i fc) = [-y/AT obs (m 7 -, M» ;fc ) + 0~75 + 
1]/Vp.ff (mi) using the approximation given by iGehreld 
(1986). The total Poisson error at each mass bin is the 
quadratic sum of all magnitude bins: 



(8) 



The errors associated with sample variance are estimated 
by accounting for the difference in the observed number 
between the two fields as: 



fe ) = \r M (M^ m ) - < / (M,, m )|/2 (9) 
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Fig. 6. — The total observed number of galaxies in UV magnitud e bins compared t o those in the literature at z ~ 4 (left) and 5 
(right). Our number counts are in excellent agreement with those by Stark ct al. (2009) except the very faint end. This discrepancy is 
likely explained by the fact that the data set (v2.0 release) use d in this work has a slightly higher completeness than theirs (vl.O release). 
Assuming the best-fit parameters of the Bouwcns ct al. (2007) determination of the UV luminosity function, we compute the expected 
galaxy number counts using Equation [3] We find good agreement between the observed and expected number counts, which confirms that 
the faint-end slope of the UV luminosity function is indeed steep. Also shown as dotted lines in each panel are the expected number counts 
when a set of shallower slopes (a = —1.6, —1.4, —1.2 from top to bottom) are assumed while the characteristic luminosity M^ 700 is fixed to 
the same value as shown in top left corner. All models are normalized at 2:350 = 24.97 and 25.71 at z ~ 4 and 5, roughly corresponding to 
the respective characteristic luminosity at the median redshift of each sample (shown as filled arrows). 



The total error budget is obtained as a quadratic sum of 
the random error and sample variance at each mass bin, 

and is given as A0 M = ^[A^] 2 + [A0f/] 2 . 

Using the prescription described above, we computed 
the stellar mass function of LBGs at z ~ 4 and 5. The 
results are shown in Figure [7] t ogether with other mea- 
surements in the li t erature (Qv iarches ini et al.ll2009l [20101: 
Stark et all l2009t [Caputi et all 120111: IGonzalez et al.l 
2011i|Saritini et al J [20 111) . We~also indicate the number 
density measured in the GOODS-S and GOODS-N fields 
as small open symbols at each mass bin (circles and trian- 
gles). Our m easurements at z y 4 are in br oad agreement 
with those of lMarchesini et al.l (|2009ll2010l blue open cir- 
cles using set 7 in their Table 4), who measured the total 
stellar mass function of all galaxies at 3 < z < 4. The 
fact that, at the massive end, their estimated number 
density is a factor of 2 — 3 higher than ours is consistent 
with their observation that only w 30% of their sample 
at the very massive end would be selected as Lyman- 
break galaxies, while the rest is either pass ive or heavily 
obscu red by dust. The measurements of iSantini et al.l 
(|2011l green downward triang les) at 3.5 < z < 4.5 , 
mostly consistent with those of lMarchesini et all {2009), 
are also in line with our data in their intermediate-mass 
bins. In their least massive bin pa 1O 1O M0, their point 
is higher by a factor of 2, which may be due to cosmic 
variance as the bin includes only a handf ul of sources. 
We al so compare our results with those of iCaputi et al.l 
inn, cyan open and filled diamonds for z ~ 3.2 and 
z ~ 3.9, filled diamonds for z ~ 5) after we corrected 
their masses by 0.16 dex. The correction is made to 
convert their masses (assuming Salpeter IMF) to ours 
(Chabrier IMF). While at log[M star /M Q ] > 10.5, their 
data are consistent with our estimates, they are system- 
atically lower by more th an a factor of 2 than our esti- 
mat es as well as those o f lMarchesini et al.l ((20091 12010D 
and ISantini et al.l (|201lD at lower masses at both red- 
shift bins. While the origin of this discrepancy is un- 
clear, possible causes include that 1) some of real high- 



rcdshift sources may have been missed by their photo- 
metric redshift selection, or 2) their incompleteness cor- 
rection may have been underestimated. Relatively large 
discrepancies present among simil ar studies based on 
the rest-frame optical selection fe.g.JFontana et al.l | 2006 ; 
Perez-Gonzalez et al.ll2008HMarchesini et al.ll2009L I201C : 



Caputi et al.l 1201 It ISantini et al.ll2011l ) may also be ex- 
plained at least in part by nonnegligible cosmic variance. 
Even the lowest-mass galaxies selected by these stud- 
ies are amongst the most massive galaxies at high red- 
shift, which ar e observed to be strongly clustered (e.g., 
ILin et alllMTI) . 

We also compare our results with those that mea- 
su red the SMF for s tar-forming galaxies. The results 
of iStark et all (|2009l . blue triangles) are in god agree- 
ment with ours when we confine our analyses only to 
galaxies at M1700 < —20 as they did. A slight discrep- 
ancy at the massive end is likely due to our stringent 
selection, which resulted in removal of a few galaxies at 
the massive end, and possibly due to slightly different 
criteria applied to remove likely low-redshift interlopers. 
These "cleaning" processes affect the massive end the 
most, and clearly, the two measurements are comparable 
in lower-mass bins at log[M s t ar /M ] < 10.5. Finally, 
IGonzalez et al.l (|2011l orange upward triangles) found 
the low-mass end slope of the SMF to be a at ~ —1.4 in 
excellent agreement with our analyses. At z ~ 3.7, their 
measurements are consistent with ours in their low-mass 
bins (10 8 M Q < M star < 1O 95 M ) even though they went 
further down to below 10 8 M Q , whereas we chose not to 
go below 10 8 7 Af Q . At z ~ 5, their overall normaliza- 
tion is about 50% higher than ours even though the es- 
timated low- mass slope is still a = —1.39 very close to 
ours «a/ = —1.36 (see discussions later). At the mas- 
sive end, however, their data points are lower by more 
th an a factor of 2 com pared to our own results and those 
of IStark et al.l (|2009D . Because their sa mple is based on 
the WFC3 Early Release Science (ERS: I Windhorst et al.l 
120111 ) data which covers about one third of the GOODS- 
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Fig. 7. — The stellar mass functions or LBGs at z ~ 4 and 5 are shown. The number density of galaxies at each stellar mass bin is 
obtained by first correcting the observed number for the selection incompleteness (for the UV selection), then summing up the contributions 
from all magnitude bins. The error bars fully account for random Poisson errors as well as those arising from sample variance. Filled circles 
represent our best estimates combining the statistics of both GOODS fields, while small open symbols ( circle and triangle ) at each mass bin 
indicate th e number density computed in each field. Also shown are sim ilar measures in the literature (Stark ct al. 2009; Marchcsini ct al. 
2009, 2010; Caputi ct al. 2011; Gonzalez ct al. 2011; Santini et al. 2011, in blue filled triangles, blue open circles, cyan diamonds, o range 
upward triangles, and green downward triangles, respectively). For z ~ 4, we show the model predictions from Guo ct al. (2011], solid 
and dotted yellow lines) based on Millennium and Millennium II simulations. The best-fit Schechter parameters (shown in solid red lines) 
indicate that the low-mass-end slope of the stellar mass function is ss —1.4, and thus is considerably shallower than the faint-end slope of 
the UV luminosity function derived for the same galaxies, —(1.7— 1.8). Top right insets illustrate the confidence intervals at the 68, 95, 
and 99% level. 



S field (i.e., ~ 15% of th e total area covered by our study 
and that of iStark et all (f2009h who used both GOODS- 
N and GOODS-S), it is possible that their massive bins 
suffer more from cosmic variance than others. 

Based on our measurements, we determined the best- 
fit parameters when a Schechter function is assumed for 
the stellar mass function: 

<MM Btar ) = (In 10)0^ [10( m -"- m -2(i+«sm F )] 

x expI-lO^"-"*-* 5 ] (10) 

Note that all parameters have subscripts "M" to distin- 
guish from similar Schechter parameters for the UV lu- 
minosity function discussed previously. To properly ac- 
count for the measurement errors in mass, we use the 
probability distribution discussed in £14.21 (also shown in 
Figure |4] and [5]) to convert the intrinsic SMF to the ob- 
served one: 



</>M(M* jm ) = Sfc0M(M» ife ) p(M^ m \M^k) 



(11) 



where M* ; fc is the true stellar mass of a galaxy and -M* jm 
is the measured one. We use the chi-square minimiza- 
tion method to search for the best parameter set (o;smf, 
M* tax ). For a given parameter set (cusmf, M* taI ), we de- 
termine the normalization parameter <jf M such that the 
total number density of galaxies at M sta r > A^min equals 
the observed number density. The minimum stellar mass 
M min is set to 10 8 - 7 M Q and 10 9 M© at z ~ 4 and 5, re- 
spectively for the reasons discussed in Sj4.2l However, our 
results are insensitive to the minimum mass as long as 
it is set within the mass range over which the determi- 
nation of stellar mass is relatively robust. On top right 
inset in Figure we show the confidence levels for the 
best-fit Schechter parameters, namely, the characteristic 
stellar mass M* tar and low-mass-end slope osmf- 

At z ~ 4, the best-fit value for the low-mass-end slope 
Qsmf is ~ —1.4. Also at z ~ 5, the best-fit slope is —1.4 



even though a large degeneracy between M* tar and ccm 
(due to lack of constraints on the massive-end) makes a 
relatively steep slope also acceptable. The slope is much 
shallower than the faint-end slope of the UV luminos- 
ity function for the same galaxies, but is similar to that 
found from stellar-mass-sele cted studies (jFontana et al.1 
120061 iMarchesini et al.ll2009f ). and a recent study of L BGs 
based on a smaller sample by iGonzalez et al.l (|2011[ ) . To 
visually illustrate this difference, we show in Figure [Jj 
(dashed lines) the best-fit Schechter function when the 
slope «smf is fixed to the faint-end slope of the UV lu- 
minosity function from lBouwens et ahl (|2007l ). the values 
of which are indicated in Figure [Jj Clearly, a very steep 
low-mass end slope «smf = —(1-7 — 1.8) does not pro- 
vide a good fit to the observed stellar mass function even 
for the z ~ 5 case. The shallow slope on the low-mass 
end seems to be a robust feature even considering some- 
what large uncertainties for the normalization parameter 
as evidenced by a relatively large cosmic variance partic- 
ularly at z ~ 5. In Table [5J we summarize the best- fit 
Schechter parameters for both samples. For the V&qq- 
band dropouts, we also provide the best-fit parameters 
when the characteristic mass M* tar is forced not to ex- 
ceed that at z ~ 4. Such a treatment generally helps 
constrain the low-mass end slope better by reducing the 
large degeneracy between the two parameters. While it 
remains to be seen if it is indeed the case, the brighten- 
ing of th e characteristic lumin osity L^y observed at high 
redshift (jBouwens et alJl2007| ) likely supports such a sce- 
nario unless the mass-to-light ratios change significantly 
with redshift to reverse the effect. We also note that re- 
moving the least robust mass bin also makes the slope 
even shallower (albeit leaving larger overall uncertain- 
ties). Finally, we c ompare our meas urements against re- 
cent predictions bv lGuo et al.l (|2011| ) using semi-analytic 
models implemented on the Millennium and Millennium 
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II simulations at z ~ 4 (shown as black lines on the 
left panel of Figure [7]) . The dotted line indicates their 
predictions for the intrinsic SMF, while the solid line rep- 
resents the values accounting for measurement errors (by 
convolving the intrinsic SMF assuming 0.3 dex scatter in 
stellar mass). On the massive end, our estimates are in 
excellent agreement with theirs. However, at the mass 
regime below the characteristic mass ~ 1O 1O M0, their 
model clearly overpredicts the abundance by a factor of 
2 — 3 compared to all the existing studies. Furthermore, 
the predicted slope for the low-mass end of the SMF is 
more in line with that of the UVLF shown in dashed line, 
much steeper than our measurements. 

Finally, we note that not all studies hnd the low-mass 
en d slopes to be as shallo w as our results as well as those 
of I Gonzalez et all (|2011f ). For example, Kaiisaw a et al.l 
(2009) studied the galaxy SMF out to z ~ 3, and re- 
ported that the low-mass end slopes become steeper 
with redshift, from ckm = —(1-2 — 1.3 ) at z ~ 0.7 to 
olm = —1.6 at z ~ 3. iSantini et al.l ()201lD reported 
a similar steepening of the slope with redshift, from 
a M = -1.44 ± 0.03 at z ~ 0.8 to -1. 86 ± 0.16 at 
z ~ 3. For galaxies at z ~ 3.2, 3.9, and 4.7. ICaputi et al.l 
( 2011) found the low -mass end slope to be « —1.8 at 
all probed redshift bins. With the current observational 
constraints, these results are not necessarily contradic- 
tory to shallower slopes reported for the SMFs of star- 
forming galaxies. Most studies based on mass-selected 
samples agree that the characteristic mass M* tar lies well 
above 10 u Mq within a range of population synthesis 
models and assumptions (see all the references discussed 
above). On the other hand, the number density of star- 
forming galaxies at around 10 n M Q is observed to be 
very low, and consequently, the SMF measured only for 
star-forming galaxies tends to have a much lower char- 
acteristic mass (10 10 — 10 10 - 6 M Q ) in most studies. Such 
low characteristic mass is found independent of whether 
the mass es are fit for individu al galaxies (such as this 
work and IGonzalez et al.l 120111 ) or by an average mass- 
to-light ratio was a ssumed (e.g... IReddv fc Steidell 120091 : 
IMcLure et~aTI 120091 ). Hypothetical^, if the fraction of 
star-forming (therefore, UV-visible) galaxies increases to- 
wards lower masses, then it is conceivable that the true 
mass function will have a high characteristic mass above 
10 ll M & but will converge towards that of SMF at low 
masses where galaxies are dominated by star-forming 
galaxies. To test this possibility, we tried to fit our data 
to a Schechter function after forcing the characteristic 
mass to be greater than 10 11 M Q . When we only in- 
clude our data points down to ps 3 x 10 9 Mq, which is 
the lowest mass bin currently probed by deep surveys 
(|Kaiisawa et al.l [2009T ) . we find that the slope as steep 
as —(1.8 — 1.9) can provide a good fit. However, once 
we include the lower mass bins, such a steep slope is no 
longer acceptable. In other words, the steep slope found 
by mass-selected studies may be explained by a combina- 
tion of the M* tar — aM degeneracy (see inset of Figure [7] 
for the elongated shape of the contours) and the limited 
mass range currently probed by mass-selected samples, if 
the low-mass end is primarily dominated by star-forming 
galaxies. Deeper data are critically needed to extend the 
mass range and the fraction of star-forming galaxies at 
different masses, to help better understand the observed 



discrepancies between different studies. 

6. DISCUSSIONS 

We have measured the SED properties of LBGs at 
z ~ 4 and 5, two of the largest photometric samples ever 
assembled consisting of 3088 and 987 galaxies, respec- 
tively. The use of TFIT, which can very robustly measure 
galaxy light within a fixed physical aperture over data 
with a wide range of PSF sizes, and a large set of galaxy 
simulations to quantify measurement errors and system- 
atics of TFIT photometry at different flux levels, have 
allowed us to make the most robust measurements to 
date of galaxy properties, such as stellar mass, age, and 
UV luminosity, for a large sample, and carry out a large 
statistical study of high-redshift star- forming galaxy pop- 
ulation as a whole. Based on these measurements, we 
concluded that, 1) there is a significant intrinsic scat- 
ter to the correlation between the observed UV lumi- 
nosity and stellar mass, 2) the number counts of galax- 
ies in the bins o f UV luminosity strongly support the 
previ ous studies (jBouwens et alJl2C)07t IReddv fc Steidell 
2009), that the UV luminosity function at these red- 
shifts increases steeply «uvlf ~ — (1-7 — 1.8) towards 
the faint end, and finally, 3) the number density of the 
same galaxies in the bins of stellar mass increases rather 
mildly q;smf ~ —(1-3 — 1.4) compared to the UV lumi- 
nosity function of the same galaxies. In what follows, we 
discuss these findings in further detail and their possible 
implications on how these UV-selected galaxies assemble 
their masses over cosmic time. 

6.1. The -Mi7oo--Mst ar Scaling Relation 

The locations of galaxies on the SFR-M star plane in- 
dicate typical star formation histories of the galaxy pop- 
ulation in question. If most galaxies have been assem- 
bling their mass smoothly over a long period of time, 
one would expect them to form a tight "main-sequence" 
on the SFR-Mstar plane as stellar mass is, to the zeroth 
order, the time integral of the past star formation his- 
tory. Indeed, a relatively tight main-sequence has been 
observed for star-fo rming galaxies in the local universe 
(jNoeske et al.ll2007D, and for massive star-forming galax- 
ies at z = 1 and 2 dElbaz et all 120071 120111: iDaddi et al.l 
I2007bt IRodiehiero et al.l 1201 ID . The observed scatter 
around the main sequence is 0.25 dex in all cases. For 
galaxy population at z > 2, however, the trend is still un- 
clear. While a similar (positive) correlati on between star 
formation rates and stellar mass exists (jPannella et al.l 
120091 : [Lee et al.ll2011l from stacking), many observations 
reported much larger scatter (a > 0.46 dex; IReddv et al.1 
2006; iMannucci et ail[200l iMagdis et al.ll2010D . 

As shown in Figure [H our results suggest that there is 
a significant scatter in the correlation between UV lumi- 
nosity (Afi7oo) and stellar mass. It is evident at both red- 
shifts that at a fixed stellar mass, there is a non-negligible 
fraction of galaxies that extend down to nearly 2 magni- 
tudes fainter than the "main-sequence" . It is difficult to 
use these observations to estimate the true extent of this 
faint tail because our samples are UV-selected, and thus 
are not complete in mass. It is worth noting that the 
UV selection misses a large number of very faint galaxies 
(M1700 > —19.5) as the completeness of the data de- 
creases at the faint end. If we correct for the missing 
galaxies, the faint-end tail is even more significant than 
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TABLE 2 

The Schechter parameters for the stellar mass function of LBGs at z ~ 4 and ~ 5. 
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a The characteristic mass logAf* tar is allowed to vary during the fit in the mass range 9 < logA'/* tar < 11.5 
b The characteristic mass logACJi is allowed to vary only below M* t = 10.35, the best-fit value for the 
£?435-band dropouts 



that observed here. Furthermore, as the UV selection 
systematically misses red galaxies with heavily obscured 
star formation, inclusion of such a population would fur- 
ther increase the scatter. 

The main limitation of interpreting the observed 
Mi7oo-M st ar scaling relation is the fact that UV light 
from thes e galaxies are systematica lly dimmed by inter- 
nal dust (jCalzetti et al.lll997l [2000h . At a fixed UV lu- 
minosity, the d istribution of UV ext i nction can vary sig - 
nificantly fe.g.. iBouwens et al.ll2009t iReddv et al.ll2010l i. 
For example, the la scatter in the color excess E(B — V) 
of 0.1 (or extinction at 1700A, A1700, of 1.0 mag) means 
that galaxies at a fixed UV luminosity can have SFRs 
different up to a factor of 6 (from — la to la). In turn, 
the intrinsic distribution of UV extinction can cause a 
significant scatter in the observed Mn 00 -M sta r scaling 
relation, even if the intrinsic relation between SFR and 
Af s tar is tight. To acco unt for this effect, we use the UV- 
slope measurem ents of iBouwens et all (|2009f ) . The main 
reason for using Bouwc ns et all (|2009f ) measurements in- 
stead of our own is their inclusion of the Hubble Ultra 
Deep Field and HUDF Parallel Fields in the analyses. 
These data not only provide higher signal-to-noise mea- 
surements of colors for sources detected in the GOODS 
data, but also allows to define a larger sample of UV- faint 
sources form which the color dis tribution is computed. 
Assuming the lCalzetti et all (|2000l ) extinction curve, the 
intrinsic scatter in the UV continuum slope /? (their Ta- 
ble 4) translates into the intrinsic scatter of 0.7 — 0.9 mag 
in the extinction at 1700 A, or ^1700- A more recent re- 
sult by [Castellano et al.l (|2011l ) also found a similar scat- 
ter. As for the observed Mi7oo-M star relation, we correct 
it for the selection completeness we estimated from our 
simulations. For example, at Af 170 o = —21 (—20, —19), 
the completeness is 89 (72, 34)%. Because our sample 
is selected based upon the UV colors and flux limit, it 
is fair to assume that at a fixed stellar mass correcting 
for the UV completeness is appropriate. In other words, 
the distribution of stellar mass at a fixed UV luminosity 
has no bearing on altering the intrinsic distribution we 
attempt to recover here. 

In Figure [51 we show the observed UV luminosity dis- 
tribution of galaxies at z ~ 4 at three different stellar 
mass bins, M star = 10 9 , IO 9 5 , 10 10 M o . Histograms rep- 
resent the distribution of absolute UV magnitudes for 



all galaxies (solid), and in individual fields (dashed and 
dotted for GOODS-S and GOODS-N, respectively). The 
smooth function in each panel approximates the emer- 
gent distribution of M1700 when there is a perfect corre- 
lation between SFR intrinsic UV luminosity, £1700 > and 
stellar mass and dust extinction is solely responsible for 
the spread in the observed UV luminosity. All functions 
have a standard deviation of a = 0.8 mag as discussed 
above, and are renormalized to match the peak of the 
observed distribution. We find that there is an excess of 
UV-faint galaxies at least in the lowest mass bin. The 
observed distribution (top left) already suggests the ex- 
cess although marginally. However, when the distribu- 
tion is corrected for the selection incompleteness, the 
existence of such excess becomes more apparent. At 
Mnoo ~ —19.5, the number of galaxies is higher by at 
least a factor of 2 when corrected for the selection effect. 
The selection completeness at this luminosity is still rel- 
atively high («50%), and thus our estimate should be 
robust. Towards the massive end, the discrepancy be- 
tween the observed distribution and that expected from 
UV extinction alone appears to lessen in both observed 
and corrected distribution. In particular, at the highest- 
mass bin, the observed distribution is in excellent agree- 
ment with the smooth function (except at the faintest 
end, where the completeness level is very low, and prob- 
ably not very robust). The implication is that the dust- 
corrected UV luminosity (directly proportional to SFR) 
and stellar mass are very tightly correlated at the mas- 
sive end, while there is a significant scatter in the relation 
at the low-mass end. 

6.2. Inferences from the Shape of the UV Luminosity 
Function and Stellar Mass Function 

Now we turn our attention to the stellar mass function 
of LBGs. In £15.11 we showed that the number counts of 
galaxies in the UV luminosity bin are consistent with a 
steep faint-end slope («uvlf ~ —(1-7 — 1.8)), in excellent 
agreement with other LF measures of high-redshift star- 
forming galaxies in th e literature (jBouwens et all [20071 : 
iReddy fcS tcidcl 2003). Interestingly, when we estimated 
the stellar mass function of the same galaxies f ij5.2l) . we 
found that the low-mass end slope of the mass function is 
relatively shallow, osmf = —(1.3 — 1.4), compared to the 
faint-end slope of the UV luminosity function (see Figure 
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Fig. 8. — The distribution of absolute UV magnitudes in stel- 
lar mass bins of M star = 10 9 , 10 9 ' 5 , 1O 1O M (from left). In upper 
panels, we show the observed distribution of UV magnitudes for 
all galaxies (solid histogram) and for galaxies in each field (dashed 
and dotted histogram) in the respective bins. The error bars in- 
dicate Poisson errors. The smooth function represents a normal 
distribution with standard deviation of a = 0.8 m ag. The scatter 
is cho sen to match the UV slope measurements by Bouwcns ct al. 
(2009), and its correlation with extinction l IMeurer et al. 1999). If 
dust extinction alone is responsible for the observed spread of the 
observed UV luminosity, then the function would closely reflect the 
observations. In lower panels, we show the same distribution but 
corrected for the selection incompleteness. 



[7]). The overall Mi 700 — M* correlation discussed in 
means that the majority of the galaxies on the low-mass 
end of the stellar mass function mainly populate the faint 
end of the UV luminosity function. The implication is 
that a substantial fraction of the galaxies on the faint end 
of the UVLF (which is very steep) have stellar mass lower 
than what the sensitivities of our data allow (i.e., M s t ar ~ 
10 8 7 M Q and 10 9 Af Q at z ~ 4 and 5, respectively; see 
EH). 

The relatively shallow low-mass end slope of the SMF 
has important implications for typical star formation his- 
tory of galaxies in our sample. If most galaxies have been 
forming stars in the past at comparable rates as when ob- 
served, then, to the zeroth order, the stellar mass func- 
tion of LBGs should have a similar shape as the UV 
LF as M st ar oc SFR. Introducing scatter to this rela- 
tion increases the overall normalization of the function, 
but does not alter the overall shape or the low-mass end 
slope as long as the scatter is roughly constant. What is 
needed to effectively map galaxies in the UV LF to that 
in the SMF is to assume luminosity-dependent stellar- 
mass-to-UV-light ratios. This may mean either the me- 
dian ratio M s tar/ii7oo is systematically lower for galaxies 
on the faint end of the UVLF than more UV-luminous 
ones or the M s t ar /^i7oo distribution has a long tail to- 
wards the low M sta r/ii7oo end. Both scenarios would 
effectively result in a large number of UV-faint galaxies 
to have stellar masses lower than the minimum stellar 
mass (« 1O 8 7 M0) that we deem reliable. 

In order to make more quantitative statements on these 
interpretations, we consider toy models in which the 
M star -Mi7oo scaling law is described by a mean and scat- 
ter. We begin by taking the UV LF, assign stellar masses 
to the UV-selected galaxies according to the scaling law, 
and thereby compute the stellar mass function of the 



same galaxies. This can be done by a simple integral: 

<?W(logM stax ) = / 0(Mi 7O o)P(log M star \M 170Q )dM 1 700 
JO 

(12) 

where M s t ar is in units of M Q and V(\og M sta , r \Mn 00 ) is 
a probability distribution of stellar mass at a fixed UV 
luminosity, which we assume to be a normal distribution. 
At the bright end (M1700 < — (20 — 21)), we fix the mean 
scaling law (where the probability distribution peaks) to 
what we observe in Figure [3] 

log Mstar = Tbright Mi 700 + 1.2 (13) 

We assume 7bri g ht = —0.415, the value that is consis- 
tent with the observed scaling rela tion at the bright end 
(jStark et alJl2009b iLee et al.ll20lU also see Figure [3]). A 
slope close to —0.4 means that the UV luminosity, £1700, 
scales linearly with stellar mass. Furthermore, we as- 
sume a constant scatter of 0.3 dex around the mean scal- 
ing law. While the scat ter has never bee n measured ro- 
bustly for these galaxies. lLee et all (|2011f) concluded that 
the M1700 — Mstar correlation must be rather tight on 
the bright end based on the stacking analysis of a large 
sample of L > L* (M1700 < —21.5) galaxies at z ~ 4. 
In any case, adding more scatter does little other than 
increasing the overall normalization of the SMF. This 
is because a larger number of galaxies on the faint end 
of the UVLF scatters into the massive end of the SMF 
than UV-luminous galaxies scatter into the low-mass end 
(similar to the "Malmquist bias"). As scatter does not 
affect the low-mass end slope, and we have no means to 
robustly measure this scatter observationally for the UV- 
faint galaxies, we set the scatter to be 0.3 dex throughout 
this discussion. 

First, we compute the SMF using the scaling law ob- 
served for the relatively UV-luminous galaxies (Equa- 
tion The results are shown in Figure M (dashed 
line) . Clearly, the low-mass end slope in this case is much 
steeper than the observed one. This is expected because 
the assumed scaling law implies that the UV luminosity 
and stellar mass scale linearly (M s t ar oc L1700), and thus 
the low-mass end slope of the SMF should have the iden- 
tical functional shape as the UVLF. Based on the fact 
that the assumed scaling law reproduces the bright end of 
the SMF, and that it is determined from the observations 
where the measurements have the least uncertainties, we 
determine that the single power-law mapping between 
the observed UV luminosity and stellar mass does not 
provide a good description of our observations. 

Next, we try a double power-law model in which the 
scaling law is fixed to Equation [13] at the bright end, 
but the slope changes to a different value at the faint 
end. This approach is a two-parameter model (the lu- 
minosity Mbroak at which the break takes place, and the 
power-law slope at the faint end 7f a i n t)- Over a wide 
range of Mb rca k and Yfaintj we computed the SMFs to find 
the best scaling law that describes the observed SMF. In 
Figure we illustrate the range of the scaling laws and 
their corresponding SMFs which return the reduced chi- 
square of \r < 1-2 (68% confidence level). The best-fit 
model is Mbrcak = —20.45 and 7f a int = —0.84. Clearly, 
all the acceptable scaling laws have the power-law slope 
7faint < —0.66, much steeper than that on the bright end 
(7bright ~ -0.4). 
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Alternatively, it is possible that the sc aling law has a 
steepe r but single slope as suggested by iGonzalez et all 
(|201 If) . In Figure HI we also show their best-fit scaling 
law (logMstar = -O.68M1700 - 4.49 with 0.5 dex scat- 
ter around the mean) and corresponding SMF predic- 
tion in dotted lines. Indeed, their scaling law is con- 
siste nt with the current observations in most luminosi- 
ties (|Stark et al.H200l ILee et al.ll20TTL filled squares and 
downward triangles, respectively) and is in good agree- 
ment with our estimation (blue swaths) assuming double 
power-laws. It also reproduces well the shallow slope of 
the galaxy SMF as measured independently by this work 
and theirs. The disagreement between the model SMF 
(dotted line) and their data (filled triangles in the right 
panel of Figure [9J can likely be explained by the fact 
that they constructed the SMF using individual mea- 
surements from their data and not scaling by the mean 
mass-to-light ratios as we did here. Their scaling law 
slightly under-predicts the galaxy space density at the 
intermediate masses (log[M s t a r/M Q ] = 9.5 — 10.5). More 
data are likely needed to improve the accuracies of mea- 
surements in the relevant mass ranges. Despite minor 
discrepancies, the two studies clearly agree in that the ob- 
served mass-to-light ratios is luminosity-dependent such 
that it is considerably lower by a factor of 2-4 for sub-L* 
galaxies compared to UV-bright galaxies. 

6.3. The Impact of Luminosity- Dependent Extinction 

As discussed in §6.21 the difference in the observed 
slopes between the UVLF and SMF of the same galax- 
ies seems robust. As shown in Figure |H1 if galax- 
ies obey a single power-law (dashed lines) on the ob- 
served M1700 — -Mgtar plane, the abundance of low-mass 
(<; 1O 9 5 M ) galaxies would be overpredicted by a fac- 
tor of several. However, the observed difference may not 
entirely be attributed to changing mass-to-light ratios if 
internal dust extinction strongly depends on UV lumi- 
nosity. If the extinction systematically increases with 
UV luminosity, the space density of UV-bright galax- 
ies would be boosted more substantially than that of 
UV- faint galaxies after correcting for dust. As a result, 
the "dust-corrected" UV luminosity function (intrinsic 
UVLF, hereafter) would have a shallower faint-end slope 
than the observed UVLF. Hence, it is important to quan- 
tify how much of this effect would be needed to explain 
the observed Mi 700 — Vf sta r scaling law without invoking 
changing mass-to-light ratios. 

Unfortunately, existing works are somewhat conflict- 
ing on this topic at pr esent ranging from a rela- 
tively strong d e pendence (IBouwens et al.l 120091 l2011at 
ILee et al.l 120111: ICastellano et al.l 120111) to little or no 
dependence (lOuchi et al.1 l2004at iDunlop etall 120111 : 
iFinkelstein et al.1 1201 lh " at z > 4. The discrepancies 
among different works likely lie in the limited knowledge 
of the precise redshifts for the samples, measurement er- 
rors arising from the use of a limited number of photo- 
metric bands (often one or two colors), and accounting 
for the selection effects of high-redshift samples that may 
mimic the trend, or quite possibly, a combination of these 
effects. Obviously, the strongest luminosity-dependent 
trend for extinction will have the most dramatic impact 
on the interpretations of our results. Therefore, we ex- 
plore this possibility by using the most recent results of 
iBouwens et al.l (|2011aj) . who reported a clear trend of 



luminosity-dependent extinction. 

We use the i r upd ated best-fit formula given by 
IBouwens et all ^0U3) . = -0.1 0(M T tv + 19.5) - 1.9 8 
and the la scatter of ap = 0.3 (|Bouwens et al.l[2"009f ). 
Assuming the correlation between the extinction (Aigoo) 
and UV slope (0) ob served for local starburst galaxies 
(jCalzetti et al.l [2000L A 1600 = 2.31/3 + 4.85), we esti- 
mate the effect of luminosity-dependent extinction on the 
shape of the UVLF. We start with a Schechter function 
for the intrinsic UVLF, then apply luminosity-dependent 
dust extinction to repopulate galaxies on the observed 
UV luminosity space, and compute the observed UVLF. 
The Schechter parameters for the intrinsic UVLF were 
varied until we reprod uced the observed UVLF at z ~ 4 
(|Bouwens et al.l 120071) . The use of a Schechter function 
is not essential to our results as long as the faint-end is 
described as a power-law. The best-fit faint-end slope 
for the intrinsic UVLF is —1.62, somewhat shallower 
than that of the observed UVLF of —1.76. Assuming 
a slightly different correlation from lMeurer et ail (1999) 
- ^1600 = 1-99/? + 4.43 - does not change the slope for 
more than Aq = 0.01. 

As expected, the resultant faint-end slope of the in- 
trinsic UVLF is shallower than that observed if there 
is a positive correlation between extinction and UV lu- 
minosity. However, it seems that the "flattening" that 
occurs as a result of luminosity-dependent extinction is 
relatively mild even for the strongest trend reported so 
far, and thus is insufficient to explain the observed shal- 
low slope -(1.3 - 1.4) for the SMF at the 1.5<t level (the 
discrepancy is at the w 2a level in the absence of such a 
trend). Hence, we conclude that changing mass-to-light 
ratios are indeed needed to explain the observed shape of 
the stellar mass function even if extinction is luminosity- 
dependent. Careful analyses on this topic will shed fur- 
ther light on more proper interpretations of our results. 

6.4. The Origin of the Changing Mass-to-Light Ratios 

In what follows, we consider mainly three physical 
scenarios that might have caused the observed scaling 
law to deviate from a single power-law. The first sce- 
nario we explore is that more UV-luminous galaxies have, 
on average, higher formation redshifts than UV-faintcr 
ones, while all galaxies have sustained their star forma- 
tion at relatively constant levels as observed. This al- 
lows us to reproduce the observed SMF while maintain- 
ing that most galaxies have formed over a long smooth 
star formation history. Such a picture is qualitatively 
in line with theoretical expectations that more massive 
galaxies formed earlier in more massive halos ("hierar- 
chical formation"). Furthermore, smooth formation his- 
tories for average galaxies are favored by numerical sim- 
ulations, as cold gas is accreted effici ently via filamen- 
tary s tructures without shock- heating (jKeres et al.ll2005L 
l2009HDekel et al.l 12009). However, this picture appears 
to directly contradict two observations. First, ages of 
galaxies would be luminosity-dependent in this scenario. 
Reading off the left panel of Figure GO galaxies with UV 
luminosity Mi 700 = —19.5 has stellar mass of roughly 
one third (blue swath) of what they would have amassed 
if the star formation varies smoothly over timd^l (dashed 

12 We note that any star formation history that varies smoothly 
over time would produce a slope of unity, including a constant SFH, 
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Fig, 9. — The inf erred M s tar--A4i700 scaling law of LBGs at z ~ 4 are sho wn t ogether with the observational measures IStark et al.l 
2009; Lee et al. 2011, cyan squares and downward triangles, respectively) (see £|6.2| , We also show the average stellar mass at a fixed UV 
luminosity directly determined from our measurements (yellow circles; see Figure [3]l. The M s tar-Mi7oo scaling law is modeled as a double 
power law, and the acceptable models that reproduce the observed stellar mass function (right ) are shown in blue. A single power-law 
model with the slope of unity well describes the observed M s tar-Mi7oo relation on the bright end (Lee et al. 2011, downward triangles), but 
produces too many low-mass galaxies as the low-mass end slope is steep, «smf = — (1-7 — 1.8), as shown in dashed lines on both panels. A 
model in which mass-to-light ratios decrease significantly towards the faint end provides a better description of th e observed number cou nts 
in stellar mass bins, such as a double power-law (blue swaths) or a steeper slope (shown in green as measured by [Gonzalez et al. 2011D, A 
steeper decline of stellar mass towards the faint end of the UVLF implies that these galaxies have not sustained their star formation at the 
rates seen at the time of observations. 



line). In other words, those galaxies should be on aver- 
age 3 times younger than those at Mi 700 = —21. In- 
stead, we find that all galaxies have similar ages ranging 
over 200 — 400 Myr when a constant star format i on his - 
tory is assumed, in agreement with iStark et al.l (|2009t) . 
However, we note that it would be hard to constrain 
population ages younger than 200 Myr as the Balmer 
break would be weak at the 0.1 — 0.2 mag level. With 
the current depth of the available data (mainly limited 
by the near-infrared data from the ground), large photo- 
metric errors can easily wash away the trend, even if it 
existed. Second, perhaps more convincing obs ervational 
constra int comes from the clustering studies. iLee et all 
(2009a|) found that, based on the measurement of corre- 
lation function of identical samples, the observed small- 
scale clustering is inconsistent with halo clustering if the 
same number density is imposed (i.e., all halos above 
a mass threshold hosts a visible galaxy). They argued 
that if only a fraction (15 — 40%) of halos hosts galax- 
ies at a given time, the shape of the correlation function 
for galaxies and halos would be well-matched. Because 
our first scenario posits that galaxies have long continu- 
ous formation histories, it directly contradicts the short 
duty cy cle implied by the clustering studies by ILee et al.l 
(|2009al) . 

The second possibility is that the formation redshifts of 
galaxies depend on luminosity, but the average star for- 
mation rates of galaxies also rise with time. Rising star 
formation history has been suggested based on several 
high-redshift observations. These include the arguments 
that galaxies at high redshift appear to be very young 
while forming a relatively tight sequence on the SFR- 

exponentially or linearly increasing/declining models unless galaxy 
ages (or formation redshifts) change with stellar mass 



M staT plane, which apparently contradict one another 
unles s rising star formation history is invok ed (|Renzinil 
120091: iMaraston et al.l 120101: ILee et al.l 120111 ) . Another 
compelling argument is that the median star formation 
rate at a fixed comoving density rises with redshift at 
least for the most UV-luminous galaxies (|Papovich et al.l 
I2011D . Rec ent numerical simul ations also support this 
view (e.g., iFinlator et all f20TTh . While a rising SFH 
would make galaxies appear young regardless of their 
formation redshifts, the second scenario also predicts a 
long duty cycle, and thus contradicts the clustering con- 
straints. 

Finally, we consider a scenario in which the average 
star formation history is more episodic. One immediate 
feature would be that galaxies are not UV-luminous at 
all times, and therefore only a fraction of galaxies would 
be observed at a given time. The exact fraction would 
be determined by how long UV-visible SF phase is sus- 
tained during the observed sampling of cosmic time. Fur- 
thermore, one would observe a substantial scatter on the 
SFR-Mgtar plane and M star --Mi7oo plane (independent of 
any variations caused by a dispersion in extinction), as 
some galaxies would be observed at the rise or decline of 
their SF episode and thus lie far outside of the "SF main 
sequence" . The episodic SFH scenario, however, does not 
preclude the existence of the main sequence albeit with 
larger scatter than that expected from continuous SF sce- 
narios. In this case, the main sequence would represent 
how the average strength of each episode depends on halo 
mass. In fact, the halo-mass-dependent UV star forma- 
tion scenario is supported by clustering studies that the 
correlation length (i.e., the clustering strength) increases 
with UV luminosity, implying that more UV-luminous 
gala xies are, on average, hosted b y more massive ha- 
los (|Giavalisco fc Dickinson! [200H lOuchi et al.l l2004bt 
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Adelberger et all 120051: TK ashika ^aet all 120061: iLee et al.l 
2006L I2009at IHildebrandt et alJ I2007L 12009ft . Such re- 
quirements can be met if more massive halos accrete, on 
average, a larger amount of cold gas to fuel star forma- 
tion episodes. In such a picture, a rising star formation 
history is a direct result of the fact that the host halo 
masses and thus possibly cold gas accretion therein in- 
crease with cosmic time. In addition, the episodic SF 
scenario would be able to explain young ages observed 
at high redshift as the galaxy light is dominated by a re- 
cently formed stellar population, even if there exist older 
generations of stars formed from earlier episodes. 

It is worth noting that, in this interpretation, the ob- 
served M s tar--^i7oo relation can be used to estimate the 
median SF duty cycle (or the duration of typical SF 
episode) . Because the unit power-law slope of the Af sta r- 
M1700 relation is achieved when SF varies smoothly over 
time, the deviation from it can immediately be inter- 
preted as due to a non-unity duty cycle. In this regard, 
the left panel of Figure [9] can provide a useful diagnos- 
tic to estimate the SF duty cycle. To the zeroth order, 
the duty cycle should be a ratio of the stellar mass that 
galaxies should have amassed if SF was continuous at the 
observed level (dashed line indicating the unit power-law 
slope) to the observed stellar mass (blue swath). Interest- 
ingly, the duty cycle of our z ~ 4 sample computed using 
this method is « 20 - 40% (evaluated at M 170 o » -19.4 
AB, roughly corresponding to the survey depth), in good 
agreement with th at estimated from our clustering study 
([Lee et al.ll2009aj ). 

One interesting implication of this interpretation is 
that the SF duty cycle changes with luminosity or halo 
mass. At the UV-luminous regime (Mi 700 < —21), star 
formation proceeds rather smoothly with long duty cy- 
cles and possibly at increasing rates. Thi s is in line with 
other studies of UV- luminous galaxies (jPapovich et al.1 
120111: ILee et al.l l201lh at these redshifts. On the other 
hand, at the UV-faint end (M1700 > —20), the dominant 
mode of SF becomes progressively shorter and thus more 
"episodic" . This may be due to more sporadic replenish- 
ment o f cold gas in relatively low-mass halos (|Lee et al.1 
l2009a| ). which likely host these UV-faint galaxies, or due 
to increasingly more efficient feedback mechanism which 
can easily shut do wn star formation by driving away 
available cold gas ()Finlator et al.ll2011f) . Whatever the 
main driver may be, the differential evolution of SF in 
dark matter halos suggested by this scenario presents a 
nontrivial test to theoretical models of galaxy formation. 

6.5. Star Formation Du ty Cycle vs. In v isible Halos 
Scenario from \Finlator et al\ \20lA ) 

In light of our findings discussed above, it is interest- 
ing to consider an alternate interpretation to the short 
SF duty cycle scenario that explains the observations at 
high redshift. The two main o bservations includ e uni- 
formly young population ages (|Stark et al.l [2009). and 
a relatively low halo occupa ncy implied by the cluster- 
ing study (|Lee et al.1 [2009a) . The former implies that 
galaxies observed at higher redshift (e.g., z ~ 6) cannot 
be progenitors of those observed at lower redshift (e.g., 
z ~ 4) because their population ages are similar. On the 
other hand, the clustering study implies that only a frac- 
tion (15 — 40%) of all halos host a visible star-forming 
galaxy at a given point in time. Both can be explained 



naturally if UV-visible star-formation does not continue 
for the Hubble time, but rather has a typical timescale 
of w 300 Myr. 

Ano ther interpretation put forward by iFinlator et al.1 
(|2011[) is that there is a significant scatter in the bary- 
onic mass and halo mass relation in low-mass halos. Mo- 
tivated by hydrodynamic simulations, they predict that 
galaxies in a significant fraction of low-mass halos will be 
affected more severely by strong outflows ( "momentum- 
driven winds" ; iDave et'a l. 2006) which suppress star for- 
mation therein. In turn, these galaxies would have much 
lower stellar mass than others hosted in halos of similar 
masses. As a result, galaxies hosted by low-mass halos 
would span a broad range of star formation rates (and 
stellar masses) , and the observed galaxies only represent 
the top end of this distribution (i.e., galaxies whose star 
formation is least suppressed by feedback) . This scenario 
can also naturally explain a low halo occupa ncy number 
similar to that measured bv lLee et al.l (|2009al 0.15—0.40) 
for LBGs, as a large fraction of halos will remain "dark" 
out of observational reach. They further predicted that 
the observed UV-selected galaxies will go on to have a 
smoothly rising star formation history with a duty cycle 
of unity, while others will remain undetected by observa- 
tions., which explains the young ages (their stellar pop- 
ulation is always dominated by recently formed stars). 

One caveat in their scenario is that the observed galax- 
ies would form a tight sequence on the SFR-Af sta r plane 
(see their Figure 10) as their SFH is smoothly rising 
with time. One direct result is that the stellar mass 
function of LBGs should still closely mirror the UV lu- 
minosity function and thus the low-mass end slope of 
the SMF is close to «uvlf ~ —2.0, a typical value 
for the faint-end sl o pe predicted by s i mulat ions (e.g., 
IFinlator et "all 120061: IChoi fc Nagaminel l2010f ). Avery 
steep SMF slope is in direct conflict with our find- 
ings and those of iGonzalez et al.l (|2011l ). that the stel- 
lar mass function rises considerably more shallowly than 
the UVLF of the same galaxies and that the stellar mass 
to UV light ratio indeed changes with luminosity. More 
careful comparison of the observed A'fstar/'Cuv with the 
predictions on a larger sample will provide a useful test 
for the validity of such a scenario. 

7. SUMMARY 

Using one of the deepest multi-wavelength data sets 
on two independent fields, we have investigated the 
statistical properties of star-forming galaxies (LBGs) at 
z ~ 4 and 5, namely, their number counts as a function 
of UV luminosity and stellar mass, as well as how these 
two quantities are related to each other. These statistics 
have important implications as to the average star 
formation history of these galaxies as they represent the 
ongoing star formation and the time integral of the past 
star formation, respectively. Based on our analyses, we 
conclude: 

1. The galaxy counts in the rest-frame UV magnitudes 
(at 1700 A) are consistent with a steep faint-end slope, 
a w —(1.7— 1.8), of the UV luminosity function, in good 
agreement with several existing studies at hig h redshift 
(jReddv fc Steidell l200l iBouwens et ail 120071 see 
and Figure [6]). 
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2. Based on the locations on the Mi7oo--Mst ar plane 
populated by galaxies in our samples, we hnd a broad 
correlation between stellar mass and UV luminosity, 
such that more UV-luminous galaxies are, on average, 
also more massive f ij6.ll see Figure [3]). However, at both 
z ~ 4 and 5, the correlation has a substantial intrinsic 
scatter, in particular, for UV-faint galaxies. This is evi- 
denced by the fact that there is a non-negligible number 
of UV-faint but massive galaxies that are present in 
our samples. Roughly 35 (20)% of the massive galaxies 
(M star > 1O 9 5 M ) are UV-faint (M uv > -20) at z ~ 4 
(5). In contrast to our findings for UV-faint galaxies, 
we also report that the region of the low stellar mass 
and high UV luminosity is largely devoid of galaxies, 
suggesting that the majority of the UV-luminous galax- 
ies may have somewhat more extended star formation 
history (for at least several hundred million years) than 
their UV- fainter counterparts. 

3. While the current data do not allow us to di- 
rectly quantify the intrinsic scatter of the SFR-Ms ta r 
scaling law, the distribution of the UV magnitudes 
at a fixed stellar mass strongly suggests that there 
exists a significant intrinsic scatter at least on the 
low-mass end (Figure [SJ. Considering the adopted 
UV color selection which systematically misses highly 
dust-obscured systems, the true distribution is likely 
even wider than that estimated in our analyses. Better 
sensitivities in the near-infrared (sampling the UV slope 
2000 — 3500A) and mid-infrared (sampling longward of 
the Balmer/4000 A break) are needed to make more pre- 
cise measurements of the intrinsic SFR-Ajf st ar scaling law. 

4. We make the statistically robust estimates of the 
stellar mass functions for LBGs at z ~ 4 and 5 ( t|5.2[l . 
Our measurements suggest that the low-mass end slope 
of the SMF is asMF ~ —(1-3 — 1.4), and thus is not 
as steep as that of the UVLF of the same galaxies at 
both reds hifts (Figure El) , in ag reement with a recent 
study by IGonzalez et al.l (|201lh based on a smaller 
sample. The direct implication is that a large fraction 
of the UV-selected galaxies are not massive enough, and 
therefore are too faint in their rest-frame optical bands, 
to be detected in the current IRAC data. The broad 
Mi7oo-M st ar correlation observed for these galaxies 
implies that those "missing" galaxies with stellar mass 
M sial ^ 10 8 ' 7 Mq are mostly UV-faint. In a scenario in 
which most galaxies have a relatively continuous star 
formation history, these galaxies would have accumu- 
lated more stellar mass (by a factor of several) than 
observed, and as a result, the SMF would more closely 
mirror that of the UVLF of the same galaxies. Hence, 
our results favor a more episodic formation history in 
which star formation rates of galaxies largely fluctuate 



over cosmic time. The duty cycle inferred from our SMF 
measurements is 20 — 40% (SJ6T4J) , in good agreement 
with th at implied from a clustering study by iLee et aLl 
(|2009al) . 

5. Using a simple toy model, we demonstrate that 
the stellar mass to UV light ratio should dec rease sys- 
tematically for galaxies with M1700 > —20.5 f £16.21 also 
see Figure [9]) in order to reproduce the observed low- 
mass end slope of the SMF. Assuming a roughly con- 
stant mass-to-light ratios for all galaxies would result in 
significant overestimation for the abundance of low-mass 
galaxies, in direct violation of the current SMF measure- 
ments. We explore a possible impact of the luminosity- 
dependent dust extinction on the interpretations of our 
results using the strongest trend reported to date, and 
conclude that it is likely insufficient to explain the ob- 
served shape of the SMF, although such a trend would 
imply a milder decline of the mass-to-light ratios towards 
faint UV luminosities than that expected in the absence 
of the luminosity-dependent extinction. We discuss sev- 
eral possible scenarios as the physical origins of the im- 
plied SFR- M s tar scaling law ( H6.4p . We conclude that 
while we cannot completely rule out a possibility that 
UV-fainter galaxies have, on average, lower formation 
redshifts based on the observations, it is most likely that 
"episodic star formation" scenario is in best agreement 
with the available set of observations including cluster- 
ing, abundance of low-mass galaxies, and population age 
constraints of high-redshift galaxies. One interesting im- 
plication of such a scenario is that the star formation 
duty cycle increases with luminosity or halo mass. A 
direct result would be that galaxies with high SFR will 
continue to form stars with duty cycles close to unity 
and assemble their mass rapidly, while galaxies with low 
SFR will lag behind as their mass assembly is punctuated 
by relatively dormant phases. The differential evolution 
of galaxies at different luminosities presents a nontrivial 
test to theoretical models of galaxy formation. 

In memory of the late Ms. Michele Dufault, a very 
bright Yale undergraduate student who would have been 
a brilliant astrophysicist of the next generation. 
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